How to use the underworld.utils function in underworld

def _calibrate_pressureField(self):
        surfaceArea = uw.utils.Integral(fn=1.0, mesh=self.mesh,
                                        integrationType='surface',
                                        surfaceIndexSet=self.topWall)
        surfacepressureFieldIntegral = uw.utils.Integral(
                fn=self.pressureField,
                mesh=self.mesh,
                integrationType='surface',
                surfaceIndexSet=self.topWall
            )
        (area,) = surfaceArea.evaluate()
        (p0,) = surfacepressureFieldIntegral.evaluate()
        offset = p0/area
        self.pressureField.data[:] -= offset
        self.pressSmoother.smooth()
        
        for material in self.materials:
            if material.viscosity:

# **Nusselt number**
# 
# The Nusselt number is the ratio between convective and conductive heat transfer
# 
# \\[
# Nu = -h \frac{ \int_0^l \partial_z T (x, z=h) dx}{ \int_0^l T (x, z=0) dx}
# \\]
# 
# 
# 
# 

# In[15]:

nuTop    = uw.utils.Integral( fn=temperatureField.fn_gradient[1], 
                              mesh=mesh, integrationType='Surface', 
                              surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])

nuBottom = uw.utils.Integral( fn=temperatureField,               
                              mesh=mesh, integrationType='Surface', 
                              surfaceIndexSet=mesh.specialSets["MinJ_VertexSet"])


# In[16]:

nu = - nuTop.evaluate()[0]/nuBottom.evaluate()[0]
print('Nusselt number = {0:.6f}'.format(nu))


# **RMS velocity**
# 

advector  = uw.systems.SwarmAdvector( swarm=swarm,       velocityField=velocityField, order=2 )


# Analysis tools
# -----
# 
# **RMS velocity**
# 
# Set up integrals used to calculate the RMS velocity.

# In[15]:

vdotv = fn.math.dot( velocityField, velocityField )
v2sum_integral  = uw.utils.Integral( mesh=mesh, fn=vdotv )
volume_integral = uw.utils.Integral( mesh=mesh, fn=1. )


# Main simulation loop
# -----
# 
# The main time stepping loop begins here. Before this the time and timestep are initialised to zero and the output statistics arrays are set up. Also the frequency of outputting basic statistics to the screen is set in steps_output.
# 
# Note that there are two ``advector.integrate`` steps, one for each swarm, that need to be done each time step.

# In[16]:

# Stepping. Initialise time and timestep.
time = 0.
step = 0
nsteps = 10
if(uw.rank()==0):

--------

        >>> mesh = uw.mesh.FeMesh_Cartesian(minCoord=(0.0,0.0), maxCoord=(1.0,2.0))
        >>> fn_1 = uw.function.misc.constant(2.0)
        >>> np.allclose( mesh.integrate( fn_1 )[0], 4 )
        True

        >>> fn_2 = uw.function.misc.constant(2.0) * (0.5, 1.0)
        >>> np.allclose( mesh.integrate( fn_2 ), [2,4] )
        True

        """
        # julian, i've dumbed this down for now as i'm not sure how to handle the
        # circular dependency.  i think the performance hit will be generally
        # negligible
        _volume_integral = uw.utils.Integral(mesh=self, fn=fn)
        return _volume_integral.evaluate()

# open hdf5 file
    h5f = h5py.File(name=filename, mode="r", driver='mpio', comm=MPI.COMM_WORLD)

    dset = h5f.get('data')
    if dset == None:
        raise RuntimeError("Can't find 'data' in file '{0}'.\n".format(filename))
    if dset.shape[1] != self.particleCoordinates.data.shape[1]:
        raise RuntimeError("Cannot load file data on current swarm. Data in file '{0}', " \
                           "has {1} components -the particlesCoords has {2} components".format(filename, dset.shape[1], self.particleCoordinates.data.shape[1]))
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    nProcs = comm.Get_size()
    
    if rank == 0 and verbose:
        bar = uw.utils._ProgressBar( start=0, end=dset.shape[0]-1, title="loading "+filename)
    
    # try and read the procCount attribute & assume that if nProcs in .h5 file
    # is equal to the current no. procs then the particles will be distributed the 
    # same across the processors. (Danger if different discretisations are used... i think)
    # else try and load the whole .h5 file.
    # we set the 'offset' & 'size' variables to achieve the above 
    
    offset = 0
    totalsize = size = dset.shape[0] # number of particles in h5 file
    
    if try_optimise:
        procCount = h5f.attrs.get('proc_offset')
        if procCount is not None and nProcs == len(procCount):
            for p_i in xrange(rank):
                offset += procCount[p_i]
            size = procCount[rank]

def solve_velocityField(self):
        fnVel = (-self._pressureField.fn_gradient + self.fn_bodyforce) * self._kMatTerm.fn

        if self._swarmVarVelocity:
            self._swarmVarVelocity.data[:] = fnVel.evaluate(self._swarm)

        if self._velocityField:
            velproj = uw.utils.MeshVariable_Projection(self._velocityField,fnVel,self._swarm)
            velproj.solve()

if not self.velocityField:
            raise ValueError("Please link a Velocity Field to the Isostasy solver")
        if not self.materialIndexField:
            raise ValueError("Please link a Material Field to the Isostasy solver")
        if not self.materialIndexField:
            raise ValueError("Please link a Material Field to the Isostasy solver")
        if not self.densityFn:
            raise ValueError("Please link a Density Field to the Isostasy solver")
        if not self.reference_mat:
            raise ValueError("Please define a reference Material for the Isostasy solver")
        
        # Create utilities
        self.MaterialIndexFieldFloat = self.swarm.add_variable(dataType="double",  count=1 )
        self.DensityVar = uw.mesh.MeshVariable(self.mesh, nodeDofCount=1)
        self.MaterialVar = uw.mesh.MeshVariable(self.mesh, nodeDofCount=1)
        self.projectorDensity = uw.utils.MeshVariable_Projection(self.DensityVar,
                                                                 self.densityFn,
                                                                 type=0 )
        self.projectorMaterial = uw.utils.MeshVariable_Projection(self.MaterialVar,
                                         self.MaterialIndexFieldFloat, type=0 )

        if not self.mesh._cself.isRegular:
            raise TypeError("You are using an irregular mesh: \
                             isostasy module only works with regular meshes")

        if self.surface is not None and not isinstance(self.surface, uw.swarm._swarm.Swarm):
            raise TypeError("'surface' must be a tuple'")
       
        self.initialized = True

"""
        # use barrier as there are some file open operations below
        # and we need to ensure that all procs have finished writing
        # before we try and open any files.
        uw.mpi.barrier()
        if uw.mpi.rank == 0:
            if not isinstance(varname, str):
                raise ValueError("'varname' must be of type str")
            if not isinstance(swarmname, str):
                raise ValueError("'swarmname' must be of type str")
            if not isinstance(filename, str):
                raise ValueError("'filename' must be of type str")
            if not isinstance(swarmSavedData, uw.utils.SavedFileData ):
                raise ValueError("'swarmSavedData' must be of type SavedFileData")
            if not isinstance(varSavedData, uw.utils.SavedFileData ):
                raise ValueError("'varSavedData' must be of type SavedFileData")
            if not isinstance(modeltime, (int,float)):
                raise ValueError("'modeltime' must be of type int or float")
            modeltime = float(modeltime)    # make modeltime a float

            # get the elementMesh - if self is a subMeshed variable get the parent
            if self.swarm != swarmSavedData.pyobj:
                raise RuntimeError("'swarmSavedData file doesn't correspond to the object's swarm")

            if not filename.lower().endswith('.xdmf'):
                filename += '.xdmf'

            # the xmf file is stored in 'string'
            # 1st write header
            string = uw.utils._xdmfheader()
            """

swarmVgrad = velocityField.fn_gradient.evaluate(swarm)
  
    stretching.data[:,0] += dt * (swarmVgrad[:,0] * stretching.data[:,0] + swarmVgrad[:,1] * stretching.data[:,1])
    stretching.data[:,1] += dt * (swarmVgrad[:,2] * stretching.data[:,0] + swarmVgrad[:,3] * stretching.data[:,1])
    
    # update plastic strain on those swarm particles that yielded
    swarmYield = viscosityFn.evaluate(swarm) < etaV
    swarmStrainRateInv = strainRate_2ndInvariantFn.evaluate(swarm)
    particleIsViscPlastic = materialVariable.evaluate(swarm) == materialV
    plasticStrainIncrement = dt * np.where(swarmYield, np.where(particleIsViscPlastic, swarmStrainRateInv, 0.0) , 0.0 )
    
    plasticStrain.data[:] += plasticStrainIncrement

    if (step%5 ==0):      
#        figVelocityPressure.save_image( outputPath + "figVP-" + str(step).zfill(4))
        projectorStress = uw.utils.MeshVariable_Projection( meshDevStress, 
                                                           fn.tensor.second_invariant(devStressFn), type=0 )
        projectorStress.solve()
#        figStrain.save_image( outputPath + "figStrain-" + str(step).zfill(4))
        
    if uw.rank()==0:   
        print('step = {0:3d}; time = {1:.3e}; xmax = {2:.3f}; pmax = {3:.4f}; cpu time = {4:.2e}'
              .format(step, time, xmax, plasticStrain.evaluate(swarm).max(), cpuTime.clock()-startTime))

    time += dt
    step += 1
  


# Post simulation images
# -----
# 

def vertical_gradient_value(self, fn=None):
        """Returns vertical gradient within the shell of scalar uw function"""

        ## Need to check here that the supplied function is
        ## valid and is a scalar. Mask function is to eliminate the
        ## effect of the values in the core.

        rGrad = uw.utils.Integral(fn * (self.unit_heightFn-self._c0)*self.maskFn, self).evaluate()[0] / self._c1

        return rGrad