How to use the underworld.mesh function in underworld

# Set size and position of dense sphere.
sphereRadius = 0.1
sphereCentre = (0., 0.7)


# Create mesh and finite element variables
# ------

# In[3]:

mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), 
                                 elementRes  = (res, res), 
                                 minCoord    = (-1., 0.), 
                                 maxCoord    = (1., 1.))

velocityField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=dim )
pressureField    = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )


# Set initial conditions and boundary conditions
# ----------
# 
# **Initial and boundary conditions**
# 
# Initialise the velocity and pressure fields to zero.

# In[4]:

velocityField.data[:] = [0.,0.]
pressureField.data[:] = 0.

import underworld as uw
import math
from underworld import function as fn

linearMesh = uw.mesh.FeMesh_Cartesian("Q2/dQ1", (32,32), (0.,0.), (1.,1.))

constantMesh = linearMesh.subMesh

velocityField = uw.mesh.MeshVariable(linearMesh,2)
velocityField.data[:] = (0.,0.)

pressureField = uw.mesh.MeshVariable(constantMesh,1)
pressureField.data[:] = 0.

# freeslip
IWalls = linearMesh.specialSets["MinI_VertexSet"] + linearMesh.specialSets["MaxI_VertexSet"]
JWalls = linearMesh.specialSets["MinJ_VertexSet"] + linearMesh.specialSets["MaxJ_VertexSet"]
freeslip = uw.conditions.DirichletCondition(velocityField, (IWalls, JWalls))


solA = fn.analytic.SolCx()

stokesSystem = uw.systems.Stokes(velocityField,pressureField,solA.fn_viscosity,solA.fn_bodyforce,conditions=[freeslip,], rtolerance=1.e-5)

# To set up Uzawa solver with mumps

def __init__(self, velocityMeshVariable, pressureMeshVariable, temperatureMeshVariable=None, pic=True, bodyforces=[], **kwargs):
        # build parent
        super(MeshVariable,self).__init__(**kwargs)
        
        # now create a bunch of stuff
        Ra=1.0
        if temperatureMeshVariable is None:
            #then build a temperature field
            mesh = velocityMeshVariable._mesh
            temperatureMeshVariable = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1, dataType="double")

        #create our own int swarm here?
        
        thermalBuoyancyTerm = ThermalBuoyancy(temperatureMeshVariable, intswarm, Ra)

        f = ForceVector(velocityMeshVariable,[])
        h = ForceVector(pressureMeshVariable,[])
        u = SolutionVector(velocityMeshVariable)
        p = SolutionVector(pressureMeshVariable)

# make directories if they don't exist
if not os.path.isdir(outputPath):
    os.makedirs(outputPath)
if not os.path.isdir(filePath):
    os.makedirs(filePath)


# Create mesh and finite element variables
# ------
# 
# Set up mesh and field variables that are solved on the mesh. See user guides for details of this process.

# In[ ]:

mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), 
                                 elementRes  = (res, res), 
                                 minCoord    = (0., 0.), 
                                 maxCoord    = (1., 1.))

velocityField       = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=2 )
pressureField       = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )
temperatureField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 )
temperatureDotField = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 )


# Set initial conditions and boundary conditions
# ----------
# 
# **Initial and boundary conditions**
# 
# Either set by perturbation function or load data from file.

def __init__(self, meshVariable, **kwargs):
        
        if not isinstance(meshVariable, mesh.MeshVariable):
            raise TypeError("'meshVariable' object passed in must be of type 'MeshVariable'")
        self._meshVariable = meshVariable
        
        # build parent
        super(Vector,self).__init__(**kwargs)

# In[ ]:

LoadFromFile = True
savedRes = 128


# **If loading from file**
# 
# Read (``savedRes`` $\times$ ``savedRes``) resolution data for $P$, $v$ and $T$ fields as well as existing summary statistics data. These are converted into lists so that the main time loop below will append with new values.
# 

# In[ ]:

if(LoadFromFile == True):
# set up mesh for savedRes*savedRes data file
    meshSaved = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), 
                                       elementRes  = (savedRes, savedRes), 
                                       minCoord    = (0., 0.), 
                                       maxCoord    = (1., 1.) )
    temperatureFieldSaved    = uw.mesh.MeshVariable( mesh=meshSaved,         nodeDofCount=1 )
    temperatureDotFieldSaved = uw.mesh.MeshVariable( mesh=meshSaved,         nodeDofCount=1 )
    pressureFieldSaved       = uw.mesh.MeshVariable( mesh=meshSaved.subMesh, nodeDofCount=1 )
    velocityFieldSaved       = uw.mesh.MeshVariable( mesh=meshSaved,         nodeDofCount=2 )

    temperatureFieldSaved.load(    inputPath+'temperatureField.h5' )
    temperatureDotFieldSaved.load( inputPath+'temperatureDotField.h5' )
    pressureFieldSaved.load(       inputPath+'pressureField.h5')
    velocityFieldSaved.load(       inputPath+'velocityField.h5')
    
    if(res==savedRes): # no remeshing needed, copy data directly
        temperatureField.data[:]    = temperatureFieldSaved.data[:]
        pressureField.data[:]       = pressureFieldSaved.data[:]

        self.specialSets["Empty"]  = lambda selfobject: uw.mesh.FeMesh_IndexSet( object           = selfobject,
                                                                                 topologicalIndex = 0,
                                                                                 size             = libUnderworld.StgDomain.Mesh_GetDomainSize( selfobject._mesh, libUnderworld.StgDomain.MT_VERTEX ))

outputPath = os.path.join(workdir,"outputs/")

if uw.rank() == 0:
    if not os.path.exists(outputPath):
        os.makedirs(outputPath)
    


# Create mesh and finite element variables
# ------
# 
# Note: the use of a pressure-sensitive rheology suggests that it is important to use a Q2/dQ1 element 

# In[5]:

mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), 
                                 elementRes  = ( resX, resY), 
                                 minCoord    = ( xmin, 0.), 
                                 maxCoord    = ( xmax, 1.),
                                 periodic    = [False, False]  ) 

velocityField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=mesh.dim )
pressureField    = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )

velocityField.data[:] = [0.,0.]
pressureField.data[:] = 0.


# ### Boundary conditions
# 
# Pure shear with moving  walls — all boundaries are zero traction with 

# 
# The mesh object has both a primary and sub mesh. "Q1/dQ0" produces a primary mesh with element type Q1 and a sub-mesh with elements type dQ0. Q1 elements have node points at the element corners, dQ0 elements have a single node at the elements centre.


mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), 
                                 elementRes  = (int(boxLength*res), res), 
                                 minCoord    = (0., 0.), 
                                 maxCoord    = (boxLength, boxHeight))

Tmesh = mesh


# Create mesh variables.  Note the pressure field uses the sub-mesh, the temperature uses the same mesh at this stage. The `temperatureDotField` variable is a work array for the advection-diffusion solver.


velocityField       = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=dim )
pressureField       = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )

temperatureField    = uw.mesh.MeshVariable( mesh=Tmesh,         nodeDofCount=1 )
temperatureDotField = uw.mesh.MeshVariable( mesh=Tmesh,         nodeDofCount=1 )

# Initialise values
velocityField.data[:]       = [0.,0.]
pressureField.data[:]       = 0.
temperatureDotField.data[:] = 0.


# Set up material parameters and functions
# -----
# 
# Set functions for viscosity, density and buoyancy force. These functions and variables only need to be defined at the beginning of the simulation, not each timestep.

nodeDofCount : int
            Number of degrees of freedom per node the variable will have

        Returns
        -------
        underworld.mesh.MeshVariable
            The newly created mesh variable.

        Example
        -------
        >>> linearMesh  = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) )
        >>> scalarFeVar = linearMesh.add_variable( nodeDofCount=1, dataType="double" )
        >>> q0field     = linearMesh.subMesh.add_variable( 1 )  # adds variable to secondary elementType discretisation
        """

        var = uw.mesh.MeshVariable(self, nodeDofCount, dataType, **kwargs)
        return var