##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
## ##
## This file forms part of the Underworld geophysics modelling application. ##
## ##
## For full license and copyright information, please refer to the LICENSE.md file ##
## located at the project root, or contact the authors. ##
## ##
##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
#
#
"""
This module contains FeMesh classes, and associated implementation.
"""
import underworld as uw
import underworld._stgermain as _stgermain
import weakref
import libUnderworld
import _specialSets_Cartesian
import underworld.function as function
import contextlib
import time
import abc
import h5py
from mpi4py import MPI
import numpy as np
[docs]class FeMesh(_stgermain.StgCompoundComponent, function.FunctionInput):
"""
The FeMesh class provides the geometry and topology of a finite
element discretised domain. The FeMesh is implicitly parallel. Some aspects
may be local or global, but this is generally handled automatically.
A number of element types are supported.
Parameters
----------
elementType : str
Element type for FeMesh. See FeMesh.elementType docstring for further info.
generator : underworld.mesh.MeshGenerator
Generator object which builds the FeMesh. See FeMesh.generator docstring for
further info.
"""
_objectsDict = { "_mesh": "FeMesh" }
_selfObjectName = "_mesh"
_supportedElementTypes = ["Q2","Q1","DQ1","DPC1","DQ0"]
def __init__(self, elementType, generator=None, **kwargs):
if not isinstance(elementType,str):
raise TypeError("'elementType' object passed in must be of type 'str'")
if elementType.upper() not in self._supportedElementTypes:
raise ValueError("'elementType' provided ({}) does not appear to be supported.\n \
Supported types are {}.".format(elementType.upper(),self._supportedElementTypes))
self._elementType = elementType.upper()
if generator == None:
if isinstance(self,MeshGenerator):
generator = self
else:
raise ValueError("No generator provided for mesh.\n \
You must provide a generator, or the mesh itself \
must be of the MeshGenerator class.")
self.generator = generator
# these lists should be populated with closure functions
# which are executed before and/or after mesh deformations
self._pre_deform_functions = []
self._post_deform_functions = []
self._arr = None
# build parent
super(FeMesh,self).__init__(**kwargs)
def _setup(self):
# add the empty set
self.specialSets["Empty"] = lambda selfobject: uw.mesh.FeMesh_IndexSet( object = selfobject,
topologicalIndex = 0,
size = libUnderworld.StgDomain.Mesh_GetDomainSize( selfobject._mesh, libUnderworld.StgDomain.MT_VERTEX ))
@property
def elementType(self):
"""
Returns
-------
str
Element type for FeMesh. Supported types are "Q2", "Q1", "dPc1" and "dQ0".
"""
return self._elementType
@property
def generator(self):
"""
Getter/Setter for the mesh MeshGenerator object.
Returns
-------
underworld.mesh.MeshGenerator
Object which builds the mesh. Note that the mesh itself may be a
generator, in which case this property will return the mesh object iself.
"""
if isinstance(self._generator, weakref.ref):
return self._generator()
else:
return self._generator
@generator.setter
def generator(self,generator):
if not isinstance(generator,MeshGenerator):
raise TypeError("'generator' object passed in must be of type 'MeshGenerator'")
if self is generator:
self._generator = weakref.ref(generator) # only keep weekref here (where appropriate) to prevent circular dependency
else:
self._generator = generator
libUnderworld.StgDomain.Mesh_SetGenerator(self._cself, generator._gen)
@property
def data_elementNodes(self):
"""
Returns
-------
numpy.ndarray
Array specifying the nodes (global node id) for a given element (local element id).
"""
uw.libUnderworld.StgDomain.Mesh_GenerateENMapVar(self._cself)
arr = uw.libUnderworld.StGermain.Variable_getAsNumpyArray(self._cself.enMapVar)
if( len(arr) % self.elementsLocal != 0 ):
raise RuntimeError("Unsupported element to node mapping for save routine"+
"\nThere doesn't appear to be elements with a consistent number of nodes")
# we ASSUME a constant number of nodes for each element
# and we reshape the arr accordingly
nodesPerElement = len(arr)/self.elementsLocal
return arr.reshape(self.elementsLocal, nodesPerElement)
@property
def data_elgId(self):
"""
Returns
-------
numpy.ndarray
Array specifying global element ids.
"""
uw.libUnderworld.StgDomain.Mesh_GenerateElGlobalIdVar(self._cself)
arr = uw.libUnderworld.StGermain.Variable_getAsNumpyArray(self._cself.eGlobalIdsVar)
return arr
@property
def data_nodegId(self):
"""
Returns
-------
numpy.ndarray
Array specifying global node ids.
"""
uw.libUnderworld.StgDomain.Mesh_GenerateNodeGlobalIdVar(self._cself)
arr = uw.libUnderworld.StGermain.Variable_getAsNumpyArray(self._cself.vGlobalIdsVar)
return arr
@property
def data(self):
"""
Numpy proxy array proxy to underlying object vertex data. Note that the
returned array is a proxy for all the *local* vertices, and it is
provided as 1d list.
As these arrays are simply proxys to the underlying memory structures,
no data copying is required.
Note that this property returns a read-only numpy array as default. If
you wish to modify mesh vertex locations, you are required to use the
deform_mesh context manager.
If you are modifying the mesh, remember to modify any submesh associated
with it accordingly.
Returns
-------
numpy.ndarray
The data proxy array.
Example
-------
>>> import underworld as uw
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) )
>>> someMesh.data.shape
(9, 2)
You can retrieve individual vertex locations
>>> someMesh.data[1]
array([ 0., -1.])
You can modify these locations directly, but take care not to tangle the mesh!
Mesh modifications must occur within the deform_mesh context manager.
>>> with someMesh.deform_mesh():
... someMesh.data[1] = [0.1,-1.1]
>>> someMesh.data[1]
array([ 0.1, -1.1])
"""
if self._arr is None:
self._arr = uw.libUnderworld.StGermain.Variable_getAsNumpyArray(self._cself.verticesVariable)
self._arr.flags.writeable = False
return self._arr
@contextlib.contextmanager
[docs] def add_post_deform_function( self, function ):
"""
Adds a function function to be executed after mesh deformation
is applied.
Parameters
----------
function : FunctionType
Python (not underworld) function to be executed. Closures should be
used where parameters are required.
"""
self._post_deform_functions.append( function )
@property
def nodesLocal(self):
"""
Returns
-------
int
Returns the number of local nodes on the mesh.
"""
return libUnderworld.StgDomain.Mesh_GetLocalSize(self._cself, 0)
@property
def nodesGlobal(self):
"""
Returns
-------
int
Returns the number of global nodes on the mesh
Example
-------
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) )
>>> someMesh.nodesGlobal
9
"""
return libUnderworld.StgDomain.Mesh_GetGlobalSize(self._cself, 0)
@property
def elementsLocal(self):
"""
Returns
-------
int
Returns the number of local elements on the mesh
Example
-------
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) )
>>> someMesh.elementsLocal
4
"""
return libUnderworld.StgDomain.Mesh_GetLocalSize(self._cself, self.dim)
@property
def elementsGlobal(self):
"""
Returns
-------
int
Returns the number of global elements on the mesh
Example
-------
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) )
>>> someMesh.elementsGlobal
4
"""
return libUnderworld.StgDomain.Mesh_GetGlobalSize(self._cself, self.dim)
[docs] def reset(self):
"""
Reset the mesh.
Templated mesh (such as the DQ0 mesh) will be reset according
to the current state of their geometryMesh.
Other mesh (such as the Q1 & Q2) will be reset to their
post-construction state.
Notes
-----
This method must be called collectively by all processes.
"""
self.generator._reset(self)
# if we have a submesh, reset it as well
if hasattr(self,"subMesh") and self.subMesh:
self.subMesh.reset()
@property
def specialSets(self):
"""
Returns
-------
dict
This dictionary stores a set of special data sets relating to mesh objects.
Example
-------
>>> import underworld as uw
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(0.,0.), maxCoord=(1.,1.) )
>>> someMesh.specialSets.keys()
['MaxI_VertexSet', 'MinI_VertexSet', 'AllWalls_VertexSet', 'MinJ_VertexSet', 'MaxJ_VertexSet', 'Empty']
>>> someMesh.specialSets["MinJ_VertexSet"]
FeMesh_IndexSet([0, 1, 2])
"""
if not hasattr(self, "_specialSets"):
class _SpecialSetsDict(dict):
"""
This special dictionary simply calls the function with the mesh object
before returning it.
"""
def __init__(self, mesh):
self._mesh = weakref.ref(mesh)
# call parents method
super(_SpecialSetsDict,self).__init__()
def __getitem__(self,index):
# get item using regular dict
item = super(_SpecialSetsDict,self).__getitem__(index)
# now call using mesh and return
return item(self._mesh())
self._specialSets = _SpecialSetsDict(self)
return self._specialSets
def _add_to_stg_dict(self,componentDictionary):
# call parents method
super(FeMesh,self)._add_to_stg_dict(componentDictionary)
componentDictionary[self._mesh.name]["elementType"] = self._elementType
def _get_iterator(self):
# lets create the full index set
iset = FeMesh_IndexSet( object = self,
topologicalIndex = 0,
size = libUnderworld.StgDomain.Mesh_GetDomainSize( self._mesh, libUnderworld.StgDomain.MT_VERTEX ) )
iset.addAll()
return iset._get_iterator()
[docs] def save( self, filename ):
"""
Save the mesh to disk
Parameters
----------
filename : string
The name of the output file.
Returns
-------
underworld.utils.SavedFileData
Data object relating to saved file. This only needs to be retained
if you wish to create XDMF files and can be ignored otherwise.
Notes
-----
This method must be called collectively by all processes.
Example
-------
First create the mesh:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) )
Save to a file (note that the 'ignoreMe' object isn't really required):
>>> ignoreMe = mesh.save("saved_mesh.h5")
Now let's try and reload. First create new mesh (note the different spatial size):
>>> clone_mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.5,1.5) )
Confirm clone mesh is different from original mesh:
>>> import numpy as np
>>> np.allclose(mesh.data,clone_mesh.data)
False
Now reload using saved file:
>>> clone_mesh.load("saved_mesh.h5")
Now check for equality:
>>> np.allclose(mesh.data,clone_mesh.data)
True
>>> # clean up:
>>> if uw.rank() == 0:
... import os;
... os.remove( "saved_mesh.h5" )
"""
if hasattr(self.generator, 'geometryMesh'):
raise RuntimeError("Cannot save this mesh as it's a subMesh. "
+ "Most likely you only need to save its geometryMesh")
if not isinstance(filename, str):
raise TypeError("'filename', must be of type 'str'")
h5f = h5py.File(name=filename, mode="w", driver='mpio', comm=MPI.COMM_WORLD)
# save attributes and simple data - MUST be parallel as driver is mpio
h5f.attrs['dimensions'] = self.dim
h5f.attrs['mesh resolution'] = self.elementRes
h5f.attrs['max'] = self.maxCoord
h5f.attrs['min'] = self.minCoord
h5f.attrs['regular'] = self._cself.isRegular
h5f.attrs['elementType'] = self.elementType
# write the vertices
globalShape = ( self.nodesGlobal, self.data.shape[1] )
dset = h5f.create_dataset("vertices",
shape=globalShape,
dtype=self.data.dtype)
local = self.nodesLocal
# write to the dset using the global node ids
dset[self.data_nodegId[0:local],:] = self.data[0:local]
# write the element node connectivity
self.data_elementNodes
globalShape = ( self.elementsGlobal, self.data_elementNodes.shape[1] )
dset = h5f.create_dataset("en_map",
shape=globalShape,
dtype=self.data_elementNodes.dtype)
if len(self.data_elgId) != len(self.data_elementNodes):
raise RuntimeError("Error in mesh.data_elementNodes - required for h5save")
local = len(self.data_elgId)
# write to the dset using the global node ids
dset[self.data_elgId[0:local],:] = self.data_elementNodes[0:local]
h5f.close()
# return our file handle
return uw.utils.SavedFileData(self, filename)
[docs] def load(self, filename ):
"""
Load the mesh from disk.
Parameters
----------
filename: str
The filename for the saved file. Relative or absolute paths may be
used, but all directories must exist.
Notes
-----
This method must be called collectively by all processes.
If the file data array is the same length as the current mesh
global size, it is assumed the file contains compatible data. Note that
this may not be the case where for example where you have saved using a
2*8 resolution mesh, then loaded using an 8*2 resolution mesh.
Provided files must be in hdf5 format, and use the correct schema.
Example
-------
Refer to example provided for 'save' method.
"""
if not isinstance(filename, str):
raise TypeError("Expected filename to be provided as a string")
# get field and mesh information
h5f = h5py.File( filename, "r", driver='mpio', comm=MPI.COMM_WORLD );
# get resolution of old mesh
res = h5f.attrs['mesh resolution']
if res is None:
raise RuntimeError("Can't read the 'mesh resolution' for the field hdf5 file,"+
" was it created correctly?")
if (res == self.elementRes).all() == False:
raise RuntimeError("Provided file mesh resolution does not appear to correspond to\n"\
"resolution of mesh object.")
dset = h5f.get('vertices')
if dset == None:
raise RuntimeError("Can't find the 'vertices' dataset in hdf5 file '{0}'".format(filename) )
dof = dset.shape[1]
if dof != self.data.shape[1]:
raise RuntimeError("Can't load hdf5 '{0}', incompatible data shape".format(filename))
if len(dset) != self.nodesGlobal:
raise RuntimeError("Provided data file appears to be for a different resolution mesh.")
with self.deform_mesh():
self.data[0:self.nodesLocal] = dset[self.data_nodegId[0:self.nodesLocal],:]
# note that deform_mesh always sets the mesh to irregular.
# reset this according to what the saved file has.
self._cself.isRegular = h5f.attrs['regular']
h5f.close()
class MeshGenerator(_stgermain.StgCompoundComponent):
"""
Abstract base class for all mesh generators.
Parameter
---------
partitioned: bool
If false, the mesh is not partitioned across entire processor pool. Instead
mesh is entirely owned by processor which generated it.
"""
_objectsDict = { "_gen": None }
_selfObjectName = "_gen"
def __init__(self, partitioned=True, **kwargs):
if not isinstance(partitioned,bool):
raise TypeError("'partitioned' parameter must be of type 'bool'.")
self._partitioned = partitioned
# build parent
super(MeshGenerator,self).__init__(**kwargs)
@property
def dim(self):
"""
Returns
-------
int
FeMesh dimensionality.
"""
return self._dim
@abc.abstractmethod
def _reset(self,mesh):
"""
Abstract class which handles mesh resetting.
"""
pass
def _add_to_stg_dict(self,componentDictionary):
# call parents method
super(MeshGenerator,self)._add_to_stg_dict(componentDictionary)
componentDictionary[self._gen.name]["partitioned"] = self._partitioned
class CartesianMeshGenerator(MeshGenerator):
"""
Abstract base class for all cartesian mesh generators.
Generators of this class provide algorithms to build meshes which are
logically and geometrically Cartesian.
Parameter
---------
elementRes: list,tuple
List or tuple of ints specifying mesh resolution. See CartesianMeshGenerator.elementRes
docstring for further information.
minCoord: list, tuple
List or tuple of floats specifying minimum mesh location. See CartesianMeshGenerator.minCoord
docstring for further information.
maxCoord: list, tuple
List or tuple of floats specifying maximum mesh location. See CartesianMeshGenerator.maxCoord
docstring for further information.
periodic: list, tuple
List or tuple of bools, specifying mesh periodicity in each direction.
"""
def __init__(self, elementRes, minCoord, maxCoord, periodic=None, **kwargs):
if not isinstance(elementRes,(list,tuple)):
raise TypeError("'elementRes' object passed in must be of type 'list' or 'tuple'")
for item in elementRes:
if not isinstance(item,(int)) or (item < 1):
raise TypeError("'elementRes' list must only contain positive integers.")
if not len(elementRes) in [2,3]:
raise ValueError("For 'elementRes', you must provide a tuple of length 2 or 3 (for respectively a 2d or 3d mesh).")
self._elementRes = elementRes
if not isinstance(minCoord,(list,tuple)):
raise TypeError("'minCoord' object passed in must be of type 'list' or 'tuple'")
for item in minCoord:
if not isinstance(item,(int,float)):
raise TypeError("'minCoord' object passed in must only contain objects of type 'int' or 'float'")
if len(minCoord) != len(elementRes):
raise ValueError("'minCoord' tuple length ({}) must be the same as that of 'elementRes' ({}).".format(len(minCoord),len(elementRes)))
self._minCoord = minCoord
if not isinstance(maxCoord,(list,tuple)):
raise TypeError("'maxCoord' object passed in must be of type 'list' or 'tuple'")
for item in maxCoord:
if not isinstance(item,(int,float)):
raise TypeError("'maxCoord' object passed in must only contain objects of type 'int' or 'float'")
if len(maxCoord) != len(elementRes):
raise ValueError("'maxCoord' tuple length ({}) must be the same as that of 'elementRes' ({}).".format(len(maxCoord),len(elementRes)))
self._maxCoord = maxCoord
self._dim = len(elementRes)
if periodic:
if not isinstance(periodic,(list,tuple)):
raise TypeError("'periodic' object passed in must be of type 'list' or 'tuple' in CartesianMeshGenerator")
for item in periodic:
if not isinstance(item,(bool)):
raise TypeError("'periodic' list must only contain booleans.")
if len(periodic) != len(elementRes):
raise ValueError("'periodic' tuple length ({}) must be the same as that of 'elementRes' ({}).".format(len(periodic),len(elementRes)))
self._periodic = periodic
for ii in range(0,self.dim):
if minCoord[ii] >= maxCoord[ii]:
raise ValueError("'minCoord[{}]' must be less than 'maxCoord[{}]'".format(ii,ii))
# build parent
super(CartesianMeshGenerator,self).__init__(**kwargs)
@property
def elementRes(self):
"""
Returns
-------
list, tuple
Element count to generate in I, J & K directions. Must be provided
as a tuple of integers.
"""
return self._elementRes
@property
def minCoord(self):
"""
Returns
-------
list, tuple
Minimum coordinate position for cartesian mesh.
Values correspond to minimums in each direction (I,J,K) of the mesh.
Note, this is the value used for initialisation, but mesh may be
advecting.
"""
return self._minCoord
@property
def periodic(self):
"""
Returns
-------
list, tuple
List of bools specifying mesh periodicity in each direction.
"""
return self._periodic
@property
def maxCoord(self):
"""
Returns
-------
list, tuple
Maximum coordinate position for cartesian mesh.
Values correspond to maximums in each direction (I,J,K) of the mesh.
Note, this is the value used for initialisation, but mesh may be
advecting.
"""
return self._maxCoord
def _add_to_stg_dict(self,componentDictionary):
# call parents method
super(CartesianMeshGenerator,self)._add_to_stg_dict(componentDictionary)
componentDictionary[self._gen.name][ "size"] = self._elementRes
componentDictionary[self._gen.name][ "minCoord"] = self._minCoord
componentDictionary[self._gen.name][ "maxCoord"] = self._maxCoord
componentDictionary[self._gen.name][ "dim"] = self._dim
if self.periodic:
componentDictionary[self._gen.name]["periodic_x"] = self._periodic[0]
componentDictionary[self._gen.name]["periodic_y"] = self._periodic[1]
if self._dim == 3:
componentDictionary[self._gen.name]["periodic_z"] = self._periodic[2]
def _reset(self,mesh):
"""
Reset method for mesh generated using the cartesian generator. This method
will reset the mesh to regular cartesian geometry.
>>> import underworld as uw
>>> import numpy as np
>>> mesh = uw.mesh.FeMesh_Cartesian(elementType='Q1')
Grab copy of original state:
>>> vertexCopy = mesh.data.copy()
Deform mesh:
>>> with mesh.deform_mesh():
... mesh.data[:] += 0.1*np.sin(mesh.data[:])
Confirm mesh is different:
>>> np.allclose(mesh.data,vertexCopy)
False
Now reset mesh and test again:
>>> mesh.reset()
>>> np.allclose(mesh.data,vertexCopy)
True
Lets do the same for the Q2 mesh:
>>> mesh = uw.mesh.FeMesh_Cartesian(elementType='Q2')
>>> vertexCopy = mesh.data.copy()
>>> with mesh.deform_mesh():
... mesh.data[:] += 0.1*np.sin(mesh.data[:])
>>> np.allclose(mesh.data,vertexCopy)
False
>>> mesh.reset()
>>> np.allclose(mesh.data,vertexCopy)
True
"""
uw.libUnderworld.StgDomain.CartesianGenerator_GenGeom( self._gen, mesh._cself, None)
mesh._cself.isRegular = True
# set algos back to regular
uw.libUnderworld.StgDomain.Mesh_SetAlgorithms( mesh._cself,
uw.libUnderworld.StgDomain.Mesh_RegularAlgorithms_New("",None) )
uw.libUnderworld.StgDomain.Mesh_DeformationUpdate( mesh._cself )
class QuadraticCartesianGenerator(CartesianMeshGenerator):
"""
This class provides the algorithms to generate a 'Q2' (ie quadratic) mesh.
"""
_objectsDict = { "_gen": "C2Generator" }
class LinearCartesianGenerator(CartesianMeshGenerator):
"""
This class provides the algorithms to generate a 'Q1' (ie linear) mesh.
"""
_objectsDict = { "_gen": "CartesianGenerator" }
class TemplatedMeshGenerator(MeshGenerator):
"""
Abstract Class. Children of this class provide algorithms to generate a
mesh by stenciling nodes on the cells of the provided geometry mesh.
Parameter:
----------
geometryMesh: underworld.mesh.FeMesh
The geometry mesh.
"""
def __init__(self, geometryMesh, **kwargs):
if not isinstance(geometryMesh,(FeMesh)):
raise TypeError("'geometryMesh' object passed in must be of type 'FeMesh'")
# note we keep a weakref to avoid circular dependencies
self._geometryMeshWeakref = weakref.ref(geometryMesh)
self._dim = self.geometryMesh.dim
super(TemplatedMeshGenerator,self).__init__(**kwargs)
@property
def geometryMesh(self):
"""
Returns
-------
underworld.mesh.FeMesh
This is the FeMesh from which the TemplatedMeshGenerator obtains
the cells to template nodes upon. Note that this class only retains
a weakref to the geometryMesh, and therefore this property may return
None.
"""
return self._geometryMeshWeakref()
def _add_to_stg_dict(self,componentDictionary):
# call parents method
super(TemplatedMeshGenerator,self)._add_to_stg_dict(componentDictionary)
componentDictionary[self._gen.name]["elementMesh"] = self.geometryMesh._cself.name
class LinearInnerGenerator(TemplatedMeshGenerator):
"""
This class provides the algorithms to generate a 'dPc1' mesh.
"""
_objectsDict = { "_gen": "Inner2DGenerator" }
def _reset(self,mesh):
"""
Reset method for mesh generated using the dPc1 generator. This method
will reset the mesh according to its geometryMesh's current state.
>>> import underworld as uw
>>> import numpy as np
>>> mesh = uw.mesh.FeMesh_Cartesian(elementType='Q2/dPc1')
Grab copy of original state:
>>> vertexCopy = mesh.subMesh.data.copy()
Deform mesh. Remember, we always modify the parent mesh:
>>> with mesh.deform_mesh():
... mesh.data[:] += 0.1*np.sin(mesh.data[:])
Confirm subMesh has been deformed:
>>> np.allclose(mesh.subMesh.data,vertexCopy)
False
Now reset mesh and test again:
>>> mesh.reset()
>>> np.allclose(mesh.subMesh.data,vertexCopy)
True
"""
uw.libUnderworld.StgFEM.Inner2DGenerator_BuildGeometry( self._gen, mesh._cself)
class dQ1Generator(TemplatedMeshGenerator):
"""
This class provides the algorithms to generate a 'dQ1' mesh.
"""
_objectsDict = { "_gen": "dQ1Generator" }
def _reset(self,mesh):
"""
Reset method for mesh generated using the dQ1 generator. This method
will reset the mesh according to its geometryMesh's current state.
>>> import underworld as uw
>>> import numpy as np
>>> mesh = uw.mesh.FeMesh_Cartesian(elementType='Q2/dQ1')
Grab copy of original state:
>>> vertexCopy = mesh.subMesh.data.copy()
Deform mesh. Remember, we always modify the parent mesh:
>>> with mesh.deform_mesh():
... mesh.data[:] += 0.1*np.sin(mesh.data[:])
Confirm subMesh has been deformed:
>>> np.allclose(mesh.subMesh.data,vertexCopy)
False
Now reset mesh and test again:
>>> mesh.reset()
>>> np.allclose(mesh.subMesh.data,vertexCopy)
True
"""
uw.libUnderworld.StgFEM.dQ1Generator_BuildGeometry( self._gen, mesh._cself)
class ConstantGenerator(TemplatedMeshGenerator):
"""
This class provides the algorithms to generate a 'dQ0' mesh.
"""
_objectsDict = { "_gen": "C0Generator" }
def _reset(self,mesh):
"""
Reset method for mesh generated using the dQ1 generator. This method
will reset the mesh according to its geometryMesh's current state.
>>> import underworld as uw
>>> import numpy as np
>>> mesh = uw.mesh.FeMesh_Cartesian(elementType='Q1/dQ0')
Grab copy of original state:
>>> vertexCopy = mesh.subMesh.data.copy()
Deform mesh. Remember, we always modify the parent mesh:
>>> with mesh.deform_mesh():
... mesh.data[:] += 0.1*np.sin(mesh.data[:])
Confirm subMesh has been deformed:
>>> np.allclose(mesh.subMesh.data,vertexCopy)
False
Now reset mesh and test again:
>>> mesh.reset()
>>> np.allclose(mesh.subMesh.data,vertexCopy)
True
"""
uw.libUnderworld.StgFEM.C0Generator_BuildGeometry( self._gen, mesh._cself)
[docs]class FeMesh_Cartesian(FeMesh, CartesianMeshGenerator):
"""
This class generates a finite element mesh which is topologically cartesian
and geometrically regular. It is possible to directly build a dual mesh by
passing a pair of element types to the constructor.
Refer to parent classes for parameters beyond those below.
Parameters
----------
elementType: str
Mesh element type. Note that this class allows the user to
(optionally) provide a pair of elementTypes for which a dual
mesh will be created.
The submesh is accessible through the 'subMesh' property. The
primary mesh itself is the object returned by this constructor.
Examples
--------
To create a linear mesh:
>>> import underworld as uw
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) )
>>> someMesh.dim
2
>>> someMesh.elementRes
(16, 16)
Alternatively, you can create a linear/constant dual mesh
>>> someDualMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) )
>>> someDualMesh.elementType
'Q1'
>>> subMesh = someDualMesh.subMesh
>>> subMesh.elementType
'DQ0'
To set / read vertex coords, use the numpy interface via the 'data' property.
"""
_meshGenerator = [ "C2Generator", "CartesianGenerator" ]
_objectsDict = { "_gen": None } # this is set programmatically in __new__
def __new__(cls, elementType="Q1/dQ0", elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.), periodic=None, partitioned=True, **kwargs):
# This class requires a custom __new__ so that we can decide which
# type of generator is required dynamically
if not isinstance(elementType,(str,tuple,list)):
raise TypeError("'elementType' object passed in must be of type 'str', 'list' or 'tuple'")
if isinstance(elementType, str):
# convert to tuple to make things easier
import re
elementType = re.split(",|-|/",elementType)
if len(elementType) > 2:
raise ValueError("A maximum of two mesh types are currently supported.")
for elType in elementType:
if not isinstance(elType,str):
raise TypeError("Items in provided 'elementType' object must be of type 'str'")
if elType.upper() not in cls._supportedElementTypes:
raise ValueError("'elementType' provided ({}) does not appear to be supported. \n \
Supported types are {}.".format(elType.upper(),cls._supportedElementTypes))
# Only supporting 2 'main' element types here
if elementType[0].upper() not in cls._supportedElementTypes[0:2]:
raise ValueError("Primary elementType must be of type '{}' or '{}'.".format(cls._supportedElementTypes[0],cls._supportedElementTypes[1]))
if len(elementType) == 2:
# all element types after first one can be in a sub-mesh.
if elementType[1].upper() not in cls._supportedElementTypes[1:]:
st = ' '.join('{}'.format(t) for t in cls._supportedElementTypes[1:])
raise ValueError("Secondary elementType must be one of type '{}'.".format(st) )
overruleDict = {}
# Only supporting 2 'main' elements types here
if elementType[0].upper() == cls._supportedElementTypes[0]:
overruleDict["_gen"] = cls._meshGenerator[0]
if elementType[0].upper() == cls._supportedElementTypes[1]:
overruleDict["_gen"] = cls._meshGenerator[1]
return super(FeMesh_Cartesian,cls).__new__(cls, objectsDictOverrule=overruleDict, **kwargs)
def __init__(self, elementType="Q1/dQ0", elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.), periodic=None, partitioned=True, **kwargs):
if isinstance(elementType, str):
# convert to tuple to make things easier
import re
elementType = re.split(",|-|/",elementType)
self._elementTypes = elementType
# ok, lets go ahead and build primary mesh (ie, self)
super(FeMesh_Cartesian,self).__init__(elementType=elementType[0], elementRes=elementRes, minCoord=minCoord, maxCoord=maxCoord, periodic=periodic, partitioned=partitioned, **kwargs)
# lets add the special sets
self.specialSets["MaxI_VertexSet"] = _specialSets_Cartesian.MaxI_VertexSet
self.specialSets["MinI_VertexSet"] = _specialSets_Cartesian.MinI_VertexSet
self.specialSets["MaxJ_VertexSet"] = _specialSets_Cartesian.MaxJ_VertexSet
self.specialSets["MinJ_VertexSet"] = _specialSets_Cartesian.MinJ_VertexSet
if(self.dim==3):
self.specialSets["MaxK_VertexSet"] = _specialSets_Cartesian.MaxK_VertexSet
self.specialSets["MinK_VertexSet"] = _specialSets_Cartesian.MinK_VertexSet
self.specialSets["AllWalls_VertexSet"] = _specialSets_Cartesian.AllWalls
def _setup(self):
# build the sub-mesh now
self._secondaryMesh = None
if len(self._elementTypes) == 2:
if self._elementTypes[1].upper() == self._supportedElementTypes[1]:
genSecondary = LinearCartesianGenerator(elementRes=elementRes, minCoord=minCoord, maxCoord=maxCoord, periodic=periodic, **kwargs)
elif self._elementTypes[1].upper() == self._supportedElementTypes[2]:
genSecondary = dQ1Generator( geometryMesh=self )
elif self._elementTypes[1].upper() == self._supportedElementTypes[3]:
genSecondary = LinearInnerGenerator( geometryMesh=self )
elif self._elementTypes[1].upper() == self._supportedElementTypes[4]:
genSecondary = ConstantGenerator( geometryMesh=self )
else:
st = ' '.join('{}'.format(t) for t in self._supportedElementTypes[1:])
raise ValueError("The secondary mesh must be of type '{}': '{}' was passed in. Tested against '{}'".format(st,self._elementTypes[1].upper(), self._supportedElementTypes[1]) )
self._secondaryMesh = FeMesh( elementType=self._elementTypes[1].upper(), generator=genSecondary )
@property
def subMesh(self):
"""
Returns
-------
underworld.mesh.FeMesh
Returns the submesh where the object is a dual mesh, or None otherwise.
"""
return self._secondaryMesh
[docs]class FeMesh_IndexSet(uw.container.ObjectifiedIndexSet, function.FunctionInput):
"""
This class ties the FeMesh instance to an index set, and stores other
metadata relevant to the set.
Parameters
----------
topologicalIndex: int
Mesh topological index for which the IndexSet relates. See
docstring for further info.
object: underworld.mesh.FeMesh
The FeMesh instance from which the IndexSet was extracted.
Example
-------
>>> amesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) )
>>> iset = uw.libUnderworld.StgDomain.RegularMeshUtils_CreateGlobalMaxISet( amesh._mesh )
>>> uw.mesh.FeMesh_IndexSet( amesh, topologicalIndex=0, size=amesh.nodesGlobal, fromObject=iset )
FeMesh_IndexSet([ 4, 9, 14, 19, 24])
"""
def __init__(self, object, topologicalIndex=None, *args, **kwargs):
if topologicalIndex not in (0,1,2,3):
raise ValueError("Topological index must be within [0,3].")
if not isinstance(object, FeMesh):
raise TypeError("'object' parameter must be of (sub)type 'FeMesh'.")
self._toplogicalIndex = topologicalIndex
super(FeMesh_IndexSet,self).__init__(object,*args,**kwargs)
@property
def topologicalIndex(self):
"""
Returns
-------
int
The topological index for the indices. The mapping is:
0 - vertex
1 - edge
2 - face
3 - volume
"""
return self._toplogicalIndex
def __repr__(self):
return repr(self.data).replace("array","FeMesh_IndexSet").replace(", dtype=uint32","")
def _checkCompatWith(self,other):
"""
Checks that these have identical topo
"""
# check parent first
super(FeMesh_IndexSet,self)._checkCompatWith(other)
if self.topologicalIndex != other.topologicalIndex:
raise TypeError("This operation is illegal. The topologicalIndex for these sets do not correspond.")
if self.object._cself != other.object._cself:
raise TypeError("This operation is illegal. The meshes associated with these IndexSets appear to be different.")
def __call__(self, *args, **kwards):
raise RuntimeError("Note that if you accessed this IndexSet via a specialSet dictionary,\n"+
"the interface has changed, and you should no longer call the object.\n"+
"This is now handled internally. Simpy use the objects directly.\n"+
"Ie, remove the '()'.")
def _get_iterator(self):
return libUnderworld.Function.MeshIndexSet(self._cself, self.object._cself)
class _FeMesh_Regional(FeMesh_Cartesian):
"""
Regional mesh class.
MinI_VertexSet / MaxI_VertexSet -> longitudinal walls : [min/max] = [west/east]
MinJ_VertexSet / MaxJ_VertexSet -> latitudinal walls : [min/max] = [south/north]
MinK_VertexSet / MaxK_VertexSet -> radial walls : [min/max] = [inner/outer]
Refer to parent classes for parameters beyond those below.
Parameter
---------
elementRes : tuple
Tuple determining number of elements (longitudinally, latitudinally, radially).
radius : tuple
Tuple determining the (inner radius, outer radius).
longExtent : float
The angular extent of the domain between great circles of longitude.
latExtent : float
The angular extent of the domain between great circles of latitude.
Example
-------
>>> (radMin, radMax) = (4.0,8.0)
>>> mesh = uw.mesh._FeMesh_Regional( elementRes=(20,20,14), radius=(radMin, radMax) )
>>> integral = uw.utils.Integral( 1.0, mesh).evaluate()[0]
>>> exact = 4/3.0*np.pi*(radMax**3 - radMin**2) / 6.0
>>> np.fabs(integral-exact)/exact < 1e-1
True
"""
def __new__(cls, **kwargs):
return super(_FeMesh_Regional,cls).__new__(cls, **kwargs)
def __init__(self, elementRes=(16,16,10), radius=(3.0,6.0), latExtent=90.0, longExtent=90.0, **kwargs):
if not isinstance( latExtent, (float,int) ):
raise TypeError("Provided 'latExtent' must be a float or integer")
self._latExtent = latExtent
if not isinstance( longExtent, (float,int) ):
raise TypeError("Provided 'longExtent' must be a float or integer")
self._longExtent = longExtent
if not isinstance( radius, (tuple,list)):
raise TypeError("Provided 'radius' must be a tuple/list of 2 floats")
if len(radius) != 2:
raise ValueError("Provided 'radius' must be a tuple/list of 2 floats")
for el in radius:
if not isinstance( el, (float,int)) :
raise TypeError("Provided 'radius' must be a tuple/list of 2 floats")
self._radius = radius
lat_half = latExtent/2.0
long_half = longExtent/2.0
# call parent cartesian mesh
# build 3D mesh centred on (0.0,0.0,0.0) - in _setup() we deform the mesh
super(_FeMesh_Regional,self).__init__(elementType="Q1/dQ0", elementRes=elementRes,
minCoord=(-long_half,-lat_half,radius[0]), maxCoord=(long_half,lat_half,radius[1]), periodic=None, **kwargs)
def _setup(self):
with self.deform_mesh():
# perform Cubed-sphere projection on coordinates
(x,y) = (np.tan(np.pi*self.data[:,0]/180.0), np.tan(np.pi*self.data[:,1]/180.0))
d = self.data[:,2] / np.sqrt( x**2 + y**2 + 1)
self.data[:,0] = d*x
self.data[:,1] = d*y
self.data[:,2] = d
#
# for index, coord in enumerate(mesh.data):
# # perform Cubed-sphere projection on coordinates
# (x,y,r) = (np.tan(np.pi*coord[0]/180.0), np.tan(np.pi*coord[1]/180.0), 1)
# d = coord[2]/np.sqrt( x**2 + y**2 + 1)
# mesh.data[index] = ( d*x, d*y, d)