How to use the underworld.function.math function in underworld

# * In steady state, if thermal energy is accurately conserved, the difference between $\langle W \rangle$ and $\langle \Phi \rangle / Ra$ must vanish, so also reported is the percentage error: 
# 
# $$ \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**

def _init_shape(self):
        center = tuple(nd(x) for x in list(self.center))
        r1 = nd(self.r1)
        r2 = nd(self.r2)
        coord = fn.input() - center
        self._fn = (fn.math.dot(coord, coord) < r2**2) & (fn.math.dot(coord, coord) > r1**2)

# a mass matrix goes into the lower right block of the stokes system coeff matrix
            self._mmatrix = sle.AssembledMatrix( self._pressureSol, self._pressureSol, rhs=self._hvector )
            # -1. as per Hughes, The Finite Element Method, 1987, Table 4.3.1, [M]

            self._compressibleTerm = sle.MatrixAssemblyTerm_NA__NB__Fn( integrationSwarm = intswarm,
                                                                        assembledObject  = self._mmatrix,
                                                                        mesh             =  mesh,
                                                                        fn               = self._fn_minus_one_on_lambda )

        if fn_stresshistory != None:
            self._NA_j__Fn_ijTerm    = sle.VectorAssemblyTerm_NA_j__Fn_ij( integrationSwarm = intswarm, 
                                                                           assembledObject  = self._fvector, 
                                                                           fn               = fn_stresshistory )
        # objects used for analysis, dictionary for organisation
        self._aObjects = dict()
        self._aObjects['vdotv_fn'] = uw.function.math.dot(self._velocityField, self._velocityField)

        super(Curvilinear_Stokes, self).__init__(**kwargs)

        for cond in self._conditions:
            if type(cond) in [uw.conditions.RotatedDirichletCondition, uw.conditions.CurvilinearDirichletCondition]:
                self.redefineVelocityDirichletBC(cond.basis_vectors)

self.bndMeshVariable.data[self.specialSets["AllWalls_VertexSet"]] = 1.0

        # note we use this condition to only capture border quadrature points
        # on the surface. For points not on the surface the bndMeshVariable will evaluate
        # <1.0, so we need to remove those from the integration as well.
        self.bnd_vec_normal  = fn.branching.conditional(
                             [ ( self.bndMeshVariable > 0.9, self._fn_unitvec_radial() ),
                               (               True, fn.misc.constant(1.0)*(1.0,0.0) ) ] )

        self.bnd_vec_tangent = fn.branching.conditional(
                             [ ( self.bndMeshVariable > 0.9, self._fn_unitvec_tangent() ),
                               (               True, fn.misc.constant(1.0)*(0.0,1.0) ) ] )

         # define solid body rotation function for the annulus

        r = fn.math.sqrt(fn.math.pow(fn.coord()[0],2.) + fn.math.pow(fn.coord()[1],2.))
        self.sbr_fn = r*self._fn_unitvec_tangent() # solid body rotation function

        self._e1 = self.add_variable(nodeDofCount=2)
        self._e2 = self.add_variable(nodeDofCount=2)

        self._e1.data[:] = self.bnd_vec_normal.evaluate(self)
        self._e2.data[:] = self.bnd_vec_tangent.evaluate(self)

        self.area = uw.utils.Integral(self.maskFn, self).evaluate()[0]
        self.full_area = 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.area
        self._c1 = uw.utils.Integral(self.maskFn*(self.unit_heightFn-self._c0)**2, self).evaluate()[0]

def _get_yieldStress3D(self):
        f = self._friction
        C = self._cohesion
        P = self.pressureField
        self.yieldStress = 6.0*C*fn.math.cos(f) + 2.0*fn.math.sin(f)*fn.misc.max(P, 0.0) 
        self.yieldStress /= (fn.math.sqrt(3.0) * (3.0 + fn.math.sin(f)))
        return self.yieldStress

# **Create variables required for plasticity calculations**

# In[ ]:

secinv = fn.tensor.second_invariant( fn.tensor.symmetric( velocityField.gradientFn ) )
coordinate = fn.coord()


# **Setup viscosity functions**
# 
# Remember to use floats everywhere when setting up functions

# In[ ]:

viscosityl1 = fn.math.exp(math.log(ETA_T)*-1.*temperatureField)
viscosityl2 = fn.math.exp((math.log(ETA_T)*-1.*temperatureField) + (1.-coordinate[1])*math.log(ETA_Y))

#Von Mises effective viscosity
viscosityp = ETA0 + YSTRESS/(secinv/math.sqrt(0.5)) #extra factor to account for underworld second invariant form

if CASE == 1:
    viscosityFn = viscosityl1
elif CASE == 2:
    viscosityFn = 2./(1./viscosityl1 + 1./viscosityp)
elif CASE == 3:
    viscosityFn = viscosityl2
else:
    viscosityFn = 2./(1./viscosityl2 + 1./viscosityp)


# **Add functions for density and buoyancy**

nbc         = negativeCond )

        self._k_term = uw.systems.sle.MatrixAssemblyTerm_NA_i__NB_i__Fn(
                        assembledObject  = K,
                        integrationSwarm = intSwarm,
                        fn = 0.5 * fn_dt * self.fn_diffusivity)

        self._m_term = uw.systems.sle.MatrixAssemblyTerm_NA__NB__Fn(
                        assembledObject  = K,
                        integrationSwarm = intSwarm,
                        fn   = 1.,
                        mesh = mesh)

        # functions used to calculate the timestep, see function get_max_dt()

        self._maxVsq  = uw.function.view.min_max(velocityField, fn_norm = uw.function.math.dot(velocityField, velocityField) )
        self._maxDiff = uw.function.view.min_max(self.fn_diffusivity)

        # Note that the c level minSep on the mesh is for the local domain
        sepFn = uw.function.misc.constant( velocityField.mesh._cself.minSep)
        minmaxSep  = uw.function.view.min_max(sepFn)
        minmaxSep.evaluate(mesh)

        self._minDx = minmaxSep.min_global()

        # the required for the solve
        self.sle = uw.utils.SolveLinearSystem(AMat=K, bVec=f, xVec=solv)

        # Check available interpolation packages

        self._mesh_interpolator_stripy = None
        self._mesh_interpolator_rbf = None