How to use the underworld.utils.Integral function in underworld

self.specialSets["surfaces_e1_normal_VertexSet"]  = lambda selfobject: surfaces_e1_normal_VertexSet

        surfaces_e2_normal_VertexSet = self.specialSets["Empty"]
        self.specialSets["surfaces_e2_normal_VertexSet"]  = lambda selfobject: surfaces_e2_normal_VertexSet

        surfaces_e3_normal_VertexSet = self.specialSets["Empty"]
        self.specialSets["surfaces_e3_normal_VertexSet"]  = lambda selfobject: surfaces_e3_normal_VertexSet

        ## mesh geometry information
        self.volume = uw.utils.Integral(self.maskFn, self).evaluate()[0]
        self.full_volume = uw.utils.Integral(1.0, self).evaluate()[0]

        ## moments of weight functions used to compute mean / radial gradients in the shell
        ## calculate this once at setup time.
        self._c0 = uw.utils.Integral(self.unit_heightFn, self).evaluate()[0] / self.volume
        self._c1 = uw.utils.Integral(self.maskFn*(self.unit_heightFn-self._c0)**2, self).evaluate()[0]


        """
        Rotation documentation.
        We create 3 basis vectors that will rotate the (x,y,z) problem to be a
        (r,lon, lat) [radial, longitudinal to colatitudinal (north)] problem.

        The rotations are performed on the local element level using the existing machinery
        provided by UW2.
        """

        # initialiase bases vectors as meshVariables

        self._e1 = self.add_variable(nodeDofCount=3)
        self._e2 = self.add_variable(nodeDofCount=3)
        self._e3 = self.add_variable(nodeDofCount=3)

def mean_value(self, fn=None):
        """Returns mean value on the shell of scalar uw function"""

        ## Need to check here that the supplied function is
        ## valid and is a scalar.

        mean_value = uw.utils.Integral(fn * self.maskFn, self).evaluate()[0] / self.volume

        return mean_value

# test it out

solver.solve( nonLinearIterate=True )
solver.print_stats()


# ## Surface integrals to calculate some average values

# In[23]:

yVelocity = velocityField[1]

surfaceArea = uw.utils.Integral(fn=1.0,mesh=mesh, integrationType='surface', surfaceIndexSet=top)
surfacePressureIntegral = uw.utils.Integral(fn=pressureField, mesh=mesh, integrationType='surface', surfaceIndexSet=top)
surfaceHeightIntegral = uw.utils.Integral(fn=yVelocity, mesh=mesh, integrationType='surface', surfaceIndexSet=top)
surfaceBaseIntegral = uw.utils.Integral(fn=yVelocity, mesh=mesh, integrationType='surface', surfaceIndexSet=base)

(area,)   = surfaceArea.evaluate()
(p0,)     = surfacePressureIntegral.evaluate() 
(vYtop,)  = surfaceHeightIntegral.evaluate()
(vYbase,) = surfaceBaseIntegral.evaluate()  

p0 /= area 
vYtop /= area
vYbase /= area 

pressureField.data[:] -= p0 



# In[24]:

v2int   = uw.utils.Integral( fn.math.dot(velocityField,velocityField), mesh )

dwint   = uw.utils.Integral( temperatureField*velocityField[1], mesh )

sinner = fn.math.dot( secinv, secinv )
vdint = uw.utils.Integral( (4.*viscosityFn*sinner), mesh )


# **Setup surface integrals used in metric functions**

# In[ ]:

rmsSurfInt = uw.utils.Integral( fn=velocityField[0]*velocityField[0], mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuTop      = uw.utils.Integral( fn=temperatureField.gradientFn[1],    mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuBottom   = uw.utils.Integral( fn=temperatureField.gradientFn[1],    mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MinJ_VertexSet"])


# **Define diagnostic functions using integrals**

# In[ ]:

def avg_temp():
    return tempint.evaluate()[0]/areaint.evaluate()[0]

def nusseltTB(temp_field, mesh):
    return -nuTop.evaluate()[0], -nuBottom.evaluate()[0]

def rms():

# 
# \\[
# 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**
# 
# The root mean squared velocity is defined by intergrating over the entire simulation domain via
# 
# \\[
# \begin{aligned}

# 
# $$ \delta = \frac{\lvert \langle W \rangle - \frac{\langle \Phi \rangle}{Ra} \rvert}{max \left(  \langle W \rangle,  \frac{\langle \Phi \rangle}{Ra}\right)} \times 100% $$

# **Setup volume integrals used in metric functions**

# In[ ]:

tempint = uw.utils.Integral( temperatureField, mesh )
areaint = uw.utils.Integral( 1.,               mesh )

v2int   = uw.utils.Integral( fn.math.dot(velocityField,velocityField), mesh )

dwint   = uw.utils.Integral( temperatureField*velocityField[1], mesh )

sinner = fn.math.dot( secinv, secinv )
vdint = uw.utils.Integral( (4.*viscosityFn*sinner), mesh )


# **Setup surface integrals used in metric functions**

# In[ ]:

rmsSurfInt = uw.utils.Integral( fn=velocityField[0]*velocityField[0], mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuTop      = uw.utils.Integral( fn=temperatureField.gradientFn[1],    mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuBottom   = uw.utils.Integral( fn=temperatureField.gradientFn[1],    mesh=mesh, integrationType='Surface', 
                          surfaceIndexSet=mesh.specialSets["MinJ_VertexSet"])


# **Define diagnostic functions using integrals**