2. Generating meshes in Florence

a) Basics of meshing and predefined shapes

In this section, we are going to have a look at generating and preparing meshes in Florence for finite element analysis. Florence offers the possibility to define meshes on a set of predefined shapes/geometries. As of version 0.1.5, Florence supports 5 types of elements. In Florence's terminology these are

1. line: Segments/beam/rod or generic 1D elements
2. tri: Triangular elements
3. quad: Quadrilateral elements
4. tet: Tetrahedral elements
5. hex: Hexahedral elements

Limited support is also available for pentagonal pent elements and generic 2D polygonal meshes. For instance, to create a mesh on a unit square domain, we simply define a Mesh instance and call the Square method on it, as follows

mesh = Mesh()
mesh.Square(nx=2, ny=2, element_type="quad")

This will create a quadrilateral mesh on a unit square domain with 4 elements. We can access the nodal coordinates and element connectivity of the mesh simply by

print(mesh.GetPoints())
print(mesh.GetElements())

Or we can simply query

print(mesh.points)
print(mesh.elements)

These attributes are standard NumPy arrays. For a more general rectangular shapes, we can define a Rectangle mesh as follows

mesh = Mesh()
mesh.Rectangle(lower_left_point=(0,0), upper_right_point=(2,1), nx=20, ny=10, element_type="tri")

This will create a triangular mesh on the specified rectangular domain with 400 elements. Finally, if we want to mesh a quadrilateral region such as parallelogram, trapezium etc, we can call the generic and rather low level map mesher from Florence called QuadrilateralProjection

mesh = Mesh().QuadrilateralProjection(c1=(0,0), c2=(2,0), c3=(2,3), c4=(1,2), npoints=10)

This will map mesh all four sides of the quadrilateral domain using 10 elements. Other predefined shapes include (not limited to)

mesh.Circle(radius=2., nrad=10, ncirc=30)

mesh.Arc(start_angle=0, end_angle=PI/2., radius=2., nrad=10, ncirc=30)

# plate with a hole
mesh.CircularPlate(side_length=30, radius=10, center=(0.,0.), ncirc=5, nrad=5, element_type="tri")

Three-dimensional meshes can be generated in a similar fashion for instance

mesh.Cube(nx=10, ny=10, nz=10, element_type="hex")

mesh.Parallelepiped(lower_left_rear_point=(0,0,0), upper_right_front_point=(2,4,10),
        nx=2, ny=4, nz=10, element_type="hex")

# for general 8 noded regions 
mesh = Mesh().HexahedralProjection(c1=(0,0,0), c2=(2,0,0), c3=(2,2,0), c4=(0,2,0.),
        c5=(0,1.8,3.), c6=(0.2,0,3.), c7=(2,0.2,3.), c8=(1.8,2,3.), npoints=6, equally_spaced=True)

mesh.Cylinder()
mesh.HollowCylinder()
mesh.Sphere()

b) Generating higher order elements

Support for arbitrary order nodal Lagrange and Lagrange Gauss-Lobatto as well as modal hierarchic elements is available. For instance, to get a 4th order Lagrange element on any type of mesh, we can call the generic function

# 4th order Lagrange Gauss-Lobatto elements
mesh.GetHighOrderMesh(p=4)
# 4th order Lagrange elements
mesh.GetHighOrderMesh(p=4, equally_spaced=True)

Apart from standard high order elements, Florence also supports generating curvilinear meshes based on the PostMesh backend. We will explore this later.

c) Mesh refinement and smoothing

Uniform mesh refinement on all kinds of meshes is supported through the generic function Refine

mesh.Refine(level=3)

Mesh smoothing (or rather non-uniform refinement) based on a given criterion is also supported. For example, to smooth all elements of a given mesh with aspect ratio 2 and above we can do

mesh.Smooth(criteria={'aspect_ratio':2})

A more sophisticated unconstrained Jacobi or Gauss-Seidal based (simple and smart) Laplacian smoothing is also available for linear and high order meshes

mesh.LaplacianSmoothing(niter=10)
mesh.LaplacianSmoothing(niter=10, algorithm="jacobi")
mesh.LaplacianSmoothing(niter=10, smart=True)
mesh.LaplacianSmoothing(niter=10, algorithm="jacobi", smart=True)

d) Mesh partitioning for parallel processing

Low and high order meshes can be partitioned for parallel processing using the built-in Partition call

partioned_mesh, partitioned_element_indices, partitioned_nodes_indices = mesh.Partition(n=12)

e) Mesh cutting operations

Low and high order meshes can be cut through using either the RemoveElements or GetLocalisedMesh functions

# geometric cuts
new_mesh = mesh.RemoveElements(xyz_min_max_array)
# algebraic cuts
new_mesh = mesh.GetLocalisedMesh(elements_indices_to_retain)

f) Mesh boolean operations

Meshes can be merged together simply through the +, - calls

mesh = Mesh().Square() + Mesh().Square(lower_left_point=(1,0))

g) Reading external meshes

The Mesh class of Florence supports reading meshes written in different formats such as Gmsh (.msh), .obj, Paraview (.vtk/.vtu), Fenics XML, .mesh, GID, Salome, Ideas (.unv), Flite (.fro), HDF5 and many more

mesh.ReadGmsh("mesh_file_name.msh", element_type="tet")
mesh.ReadOBJ("mesh_file_name.obj", element_type="tri")
# or simply call
mesh.Read(filename, reader_type="gmsh", element_type="tet")

h) Writing meshes to external formats

The Mesh class of Florence supports writing meshes to different formats such as Gmsh (.msh), .obj, Paraview (.vtk/.vtu), Fenics XML, .mesh, GID, Salome, Ideas (.unv), Flite (.fro), HDF5 and many more

mesh.WriteGmsh("mesh_file_name.msh")
mesh.WriteOBJ("mesh_file_name.obj")
mesh.WriteVTK("mesh_file_name.vtu")
mesh.WriteHDF5("mesh_file_name")
# or simply call
mesh.Write(filename, reader_type="gmsh")

i) Mesh properties and attributes

Inquiring mesh edges/faces, boundary edges/faces is as simple as

all_edges = mesh.GetEdges()
boundary_edges = mesh.GetBoundaryEdges()
interior_edges = mesh.GetInteriorEdges()

all_faces = mesh.GetFaces()
boundary_faces = mesh.GetBoundaryFaces()
interior_faces = mesh.GetInteriorFaces()

Associativity of elements with edges/faces can also be enquired using

edge_to_element_map = mesh.GetElementsWithBoundaryEdges()
face_to_element_map = mesh.GetElementsWithBoundaryFaces()

Finding which nodes are shared between which elements can be done as follows

elements_sharing_nodes = mesh.GetNodeCommonality()

A series of other attributes are available to query, such as

mesh.GetNumberOfElements()
mesh.GetNumberOfNodes()
mesh.InferElementType()
mesh.InferBoundaryElementType()
mesh.InferPolynomialDegree()
mesh.InferSpatialDimension()
mesh.InferNumberOfNodesPerElement()
mesh.InferNumberOfNodesPerLinearElement()
mesh.CreateDummyLowerDimensionalMesh()
mesh.CreateDummyUpperDimensionalMesh()
mesh.IsHighOrder
mesh.IsCurvilinear
mesh.IsEquallySpaced
mesh.IsSimilar(other_mesh)
mesh.IsEqualOrder(other_mesh)
mesh.Lengths()
mesh.Areas()
mesh.Volumes()
mesh.AspectRatios()

j) Mesh conversion routines

The Mesh class of Florence supports converting an already loaded mesh to a different type of mesh, for instance

mesh.ConvertTrisToQuads()
mesh.ConvertQuadsToTris()
mesh.ConvertTetsToHexes()
mesh.ConvertHexesToTets()

Furthermore higher order elements can also be tessellated locally to generate low order mesh with the same number of nodes in the computational mesh

mesh.ConvertToLinearMesh()