Element Mesh Generation
| Introduction | Markers |
| Approximation of Regions with ElementMesh | Element Meshes in Other Functions |
| Element Meshes with Subregions |
In order to use mesh generation functionality, the finite element method (FEM) package needs to be loaded.
Introduction
Many numerical solution techniques work by replacing a region of interest with an approximation of that region. This approximation is called a discrete region. The discrete region is partitioned into a collection of smaller elements that, as a sum, make up the entire discrete region. This partitioned discrete region is called a mesh. Finding the numerical solution is then based on computing the solution on the smaller elements and then combining the partial solutions into a solution over the entire mesh.
NDSolve, for example, internally converts the region 
into an ElementMesh object. This ElementMesh is a discrete, approximate version of the region 
over which the numerical analysis is conducted. Numerical functions such as NDSolve are capable of receiving an ElementMesh as input instead of the symbolic region description. This gives great flexibility over the mesh generation process; the element mesh, in fact, could have been generated by an external tool, for example.
There are various ways to create an ElementMesh, and various functions are provided to assist during the mesh creation process. The main function to create an ElementMesh is called ToElementMesh. ToElementMesh allows for various conceptually different methods to create an ElementMesh:
Implicit element mesh creation is based on converting an implicit function such as an ImplicitRegion to an ElementMesh. By contrast, explicit mesh generation is based on converting an explicit representation such as a GraphicsComplex into an element mesh. Manual mesh generation also falls into this explicit category. Here an explicit set of mesh elements is given to form an element mesh.
Both implicit and explicit element mesh generation can be broken up into yet smaller pieces. The function ToBoundaryMesh generates a boundary representation of an implicit or explicit input. This boundary representation can then in turn be given to ToElementMesh to form a full mesh. Among other uses, ToBoundaryMesh is useful if a numerical method only needs a boundary representation.
One important point to keep in mind is that regardless of which way the element mesh is created, the element mesh is, except in simple cases, only an approximation to the exact region. The exactness with which the element mesh captures essential parts of the region is crucial to the overall quality of the numerical solution.
Once a mesh is created, it can then be passed to numerical functions such as NDSolve.
Passing an ElementMesh to NDSolve
As a first example, a Poisson partial differential equation (PDE) with Dirichlet boundary conditions is solved in a standard way. During this process, a mesh is internally generated. In the next step, that mesh is predefined and given as an argument to NDSolve.
Note that the Interpolation function stores an ElementMesh if NDSolve made use of the finite element method. If an ElementMesh is stored with the interpolation function, then that can be extracted.
Instead of specifying an implicit parametric region, it is also possible to specify an explicit ElementMesh. This can be done by using ToElementMesh.
More information about the visualization of element meshes can be found in the element mesh visualization tutorial.
Next, the same PDE is solved, this time with only the explicit mesh defined.
One should keep in mind that four elements are not sufficient in most cases to represent an accurate solution.
Passing Options for the ElementMesh Creation to NDSolve via MeshOptions
All options for ToElementMesh and ToBoundaryMesh can be given to NDSolve directly.
As an alternative, the element mesh can be generated prior to the simulation and given to NDSolve.
More details about specifying options for NDSolve and its solution stages can be found in the Details and Options sections of NDSolve.
If a mesh is given to NDSolve, then the mesh options have no effect.
Comparing ElementMesh and MeshRegion
Before a deeper discussion of the ElementMesh functionality is given, it is instructive to compare the ElementMesh object to the MeshRegion object.
The default "MeshOrder" for an ElementMesh is 2.
The ability to deal with higher-order elements, expressed through the mesh order, has two distinct advantages:
- if needed, a conversion to e.g MeshRegion is easily possible
- simple; only one converter ToBoundaryMesh and one converter ToElementMesh
The ElementMesh is able to approximate the Disk better because it uses a second-order approximation.
The conversion between MeshRegion and ElementMesh is easy.
Using a converted MeshRegion may not give expected results when the originating ElementMesh was of higher than first order, which is the default. One way to avoid that is to convert the ElementMesh to a first-order mesh and then convert to a MeshRegion.
In order to provide the advantages compared to the MeshRegion data structure, a few caveats are to be considered. In an ElementMesh, the difference between a boundary mesh and a full mesh is indicated by the presence or absence of full mesh elements. A boundary element mesh has Automatic set for the full mesh elements.
A boundary ElementMesh does not need to be a closed curve.
This is needed for numerical algorithms (e.g. boundary integration) that do not need a full mesh. A not-closed boundary curve can, however, not be converted to a full element mesh. Note how that contrasts to a BoundaryMeshRegion: a BoundaryMeshRegion is always a representation of a full region by using the boundary to enclose that region.
A further difference between an ElementMesh and a MeshRegion is that a MeshRegion may contain lower-dimensional components that are detached from the full-dimensional mesh region elements.
Note that the lower-dimensional component was ignored during the conversion.
A boundary element mesh may have internal structure; for example, to represent two material regions.
Note how the internal structure is still present in the final ElementMesh.
Following, for completeness, is an itinerary of the advantages and disadvantages to directly using an ElementMesh. Some of the remaining items will be further discussed in this tutorial.
Approximation of Regions with ElementMesh
The methods to create an ElementMesh are by using ToElementMesh for a full mesh or ToBoundaryMesh for a boundary mesh representation.
When ToElementMesh is called, ToBoundaryMesh is first called internally. While everything can be done with ToElementMesh, it can be convenient to inspect a boundary mesh before a full mesh is generated. ToBoundaryMesh is particularly useful when working with markers on the boundary.
Manual Mesh Creation
A minimal ElementMesh is made up from coordinates and elements. The following elements are available:
- 1D: LineElement
An element is specified by its type, for example TriangleElement and a list of integer lists.
The integers are indices (also called incidents) that correspond to coordinates. In the preceding case, the TriangleElement contains two triangle elements, one made up of the incidents {1,2,3} and one made up of the incidents {3,4,1}.
To construct a mesh, coordinates must be given. The coordinates can be 1D, 2D or 3D but must match the element type. For a triangle mesh, 2D coordinates must be given.
In this example, the incidents of the second triangle element are {3,4,1} and refer to coordinates {1.,1.}, {1.,0.} and {0.,0.}, respectively.
The concept of the element incidents is closely related to that of a GraphicsComplex.
Line Meshes
For a 1D mesh, the mesh elements ei are LineElement. The boundary elements bi are then PointElement.
Triangle Meshes
For a 2D mesh, the mesh elements ei are TriangleElement or QuadElement. The boundary elements bi are then LineElement.
To create a triangle mesh, coordinates and triangle elements are needed. A linear triangle has three incidents and the element type is TriangleElement. The number of integers in an incident list determines the order of the mesh. In the triangle case, three incidents per element corresponds to linear triangles, which correspond to a first-order mesh. The incidents must be given in a counterclockwise manner.
The element-element connectivity provides neighboring information about elements.
The first edge of the first element connects nodes {3,2} and does not have a neighbor element, and is thus 0. The second edge connecting nodes 2 to the center node connects to element 2. The third element has an edge connecting the center node to node 3; this edge is connected to element 6. All subsequent elements are treated in the same fashion.
The element containers, like TriangleElement, can also hold markers. These are convenient for marking different material domains. The number of markers must be the same as the number of elements.
Until now, the triangle meshes were of first order. A quadratic triangle mesh has six incidents per element. The additional coordinates are the mid-side node. Thus, a quadratic triangle element has six incidents. The first three are the linear incidents, and the following three are the second-order incidents. See TriangleElement.
More information about the visualization of element meshes can be found in the element mesh visualization tutorial.
Quad Meshes
QuadElement meshes behave exactly the same as TriangleElement meshes, with the exception that for linear quad elements, four incidents per element are needed, and for quadratic elements, eight incidents per element are needed.
Markers can be given in exactly the same manner as in TriangleElement.
Mixed Element Type Meshes in 2D
For a 2D mesh, the mesh elements ei can be a combination of TriangleElement and QuadElement. The boundary elements bi are then LineElement.
All elements in a mixed-element mesh need to be of the same order; it is not possible to have first-order triangle elements and second-order quad elements in the same mesh.
Mixed-type element meshes can also hold markers.
The element connectivity keeps the information about marker boundaries. Whenever the marker value of two connecting elements is different, the element connectivity entry for that element is negative.
In the preceding graphic, the element number 2 with a marker 1 is connected to element number 4 with a marker 2. The element connectivity of this jump in marker values is stored in the negative sign.
An example of a mixed-element 2D mesh, which can also be used as a boundary layer mesh, can be found in the PDE model collection.
Boundary Meshes in 2D
Boundary meshes are useful to generate full meshes. In 2D, the boundary elements bi are the LineElement.
For a boundary mesh, you can add arbitrary points inside the bounding region to be part of the mesh.
Tetrahedron Meshes
Three-dimensional manual mesh creation follows the same ideas as in one and two dimensions.
Hexahedron Meshes
Mixed Element Type Meshes in 3D
For a 3D mesh, the mesh elements ei can be either a TetrahedronElement or HexahedronElement but not both.
Boundary Meshes in 3D
Boundary meshes are useful to generate full meshes. In 3D, the boundary elements bi are TriangleElement or QuadElement.
Numerical Regions
The purpose of a NumericalRegion is to combine a symbolic region representation and boundary and full element mesh. This provides additional flexibility beyond the currently automated mesh generation. Here is an example: even though there is an option to control the maximum boundary cell size via the "MaxBoundaryCellMeasure" option for ToBoundaryMesh, it is not possible to refine only parts of the boundary. A refinement of the entire boundary, however, may lead to an unnecessary large number of mesh elements. One way to work around that is to manually generate the boundary mesh. Furthermore, if a symbolic description for the region is present, a high-quality second-order mesh with minimal elements can be generated.
In this case, the mesh generation process failed to refine the small hole in the region center. While increasing the PrecisionGoal can help, it is not necessarily the best solution.
The number of mesh elements has increased significantly. One way to circumvent that is to combine the symbolic region description, manually generate a boundary mesh, and combine them in a NumericalRegion for final meshing. NumericalRegion is generated with ToNumericalRegion.
To manually generate the boundary mesh, the idea is to generate each part of the region separately and then combine them in a new boundary mesh.
The new boundary mesh is generated by joining the coordinates from both meshes and then constructing a new boundary elements list. The new boundary elements list is the combination of the boundary elements from one boundary mesh and an index-shifted version of the boundary elements from the second boundary mesh. The index shifting is necessary because the coordinates from the second mesh are now behind the coordinates from the first mesh.
Note that there are far fewer mesh elements for the same approximation quality.
Merging boundary element meshes works best if the boundary element meshes are disjoint. Intersections are typically also fine, but overlapping edges in 2D or overlapping faces in 3D can pose problems. This approach should be used if the difference of scale of the subregions is large. For other cases, the creation of subregions is better handled by Boolean operations in conjunction with the "RegionHoles" and "RegionMarker" options described in the next section.
Region Approximation Quality
For graphics primitives like Line or Polygon or a MeshRegion, a conversion to an ElementMesh is lossless. For example, an ElementMesh representation of a Rectangle is as exact or inexact, as the Rectangle itself represents a region. One way to estimate the quality of an approximation is to compare the area of the region in question, when possible.
Compare this to a Disk, for example. No matter how finely the mesh is made and how high of a mesh order the elements are, the discretization will only be an approximation.
The discrepancy is there, despite that the boundary points are close on the boundary.
In the case of a Disk, the mesh order plays a crucial role, since the ElementMesh can hold curved boundary mesh elements.
Also, the granularity of the boundary plays a role in how well a region can be approximated. The overall accuracy of the region conversion is controlled via the AccuracyGoal.
Internally, a NumericalRegion is created first. ToBoundaryMesh is called on that numerical region, and then the full mesh is created. This may introduce new nodes on the boundary. For a second-order mesh, the mid-side nodes are additionally inserted. Then a function to improve the position of these new boundary nodes is called.
In this case, the result is not different from a first-order approximation.
A NumericalRegion can hold an exact symbolic representation of the region, a boundary ElementMesh and a full ElementMesh. The improvement of the position of the boundary nodes can happen because the symbolic region is available and the exact position of the boundary nodes and also the higher-order nodes on the boundary can be found. This process is explained in more detail in the section Numerical Regions.
For the creation of the boundary mesh, two boundary mesh generators are available. The first is based on RegionPlot and called "RegionPlot". This provides a fast boundary approximation. The default boundary mesh generator is called "Continuation" and based on a continuation method. This boundary mesh generator is somewhat slower, but high accuracy can be achieved.
As a contrast, compare the result to the fast region plot method.
The visual inspection already shows that the cusps are not well resolved. The "RegionPlot" approximation may be improved by increasing the number of sample points (which is similar to PlotPoints).
Element Mesh Quality
Regardless of how well a region is approximated by an element mesh, the elements in the mesh themselves also exhibit an influence on the solution of a numerical task. For different numerical applications, different constraints on the element mesh come into play; in the finite element method, the following induce large discretization errors [4]:
The overall quality of a mesh can be expressed as a number between 
and 
, with 
being the best quality under a given quality estimate; 
indicates the worst quality. A negative quality indicates that the incidents of an element are not in the proper order, that the incidents produce a self-intersecting element, or both.
An ElementMesh has a notion of its quality. Once an ElementMesh is created, the quality of the mesh can be queried.
The result of a "Quality" computation is a quality estimate for each mesh element grouped by mesh elements.
Second-order element meshes are transformed to first-order element meshes, and then the quality is computed.
A discussion of the specific quality estimates used for different mesh elements follows.
Line Element Mesh
For each line element, the mesh quality defaults to 1.
Triangle Element Mesh
The quality of a triangle element mesh is computed according to the following formula [1, 2]:
Here 
is the area of the triangle and 
is the 
th edge length.
Quad Element Mesh
The quality of a quad element mesh is computed according to the following formula:
Here 
is the area of the quad and 
is the 
th edge length.
Tetrahedron Element Mesh
The quality of a tetrahedron element mesh is computed according to the following formula [1, 3]:
Here 
is the volume of the tetrahedron and 
is implemented as 
, where 
is the 
th edge length.
Hexahedron Element Mesh
The quality of a hexahedron element mesh is computed according to the following formula:
Here 
is the volume of the hexahedron and 
is the 
th edge length.
References
[1] J. R. Shewchuk, "What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measures (Preprint)," 2002, unpublished.
[2] R. P. Bhatia and K. L. Lawrence, "Two-Dimensional Finite Element Mesh Generation Based on Stripwise Automatic Triangulation," Computers and Structures, 36, 1990 pp. 309–319.
[3] V. N. Parthasarathy, C. M. Graichen, and A. Hathaway, "Fast Evaluation & Improvement of Tetrahedral 3-D Grid Quality," 1991 (manuscript).
[4] S.-W. Cheng, T. K. Dey, and J. R. Shewchuk, Delaunay Mesh Generation, Boca Raton: CRC Press, 2013.
Visualize Low-Quality Elements
When a mesh is generated that has low-quality elements, it may become necessary to visualize those elements and then address the quality in that specific area.
While the average quality of a mesh can be improved, it is possible that elements with a bad quality estimate are in corners of the geometry. Those cannot always be improved because the geometry dictates the shape. Nevertheless, it is important to be aware of them.
It can be noted that the overall average quality has improved. The mesh now, however, has more elements, which will result in a slightly longer computation time for a numerical analysis.
To visualize the elements that are below a certain threshold, the poor-quality elements are selected from the mesh element quality list and visualized.
Element Meshes with Subregions
It is common for a PDE to interact with a region that is made up of multiple materials. The solution of PDEs will be of a higher quality if the mesh elements do not cross the internal boundaries. To illustrate this, a PDE with a variable diffusion coefficient is reconsidered and solved over two regions. (See: Solving Partial Differential Equations with the Finite Element Method.) One region respects the internal boundary, while the other does not.
The diffusion coefficient has a jump discontinuity at 
.
The next example shows a circular region with a subregion and holes inside the subregion.
To easily create a mesh from the region that respects the interior boundary, the region is set up in such a way that all interior boundaries are created.
In the next step, what is a region hole and what is not is inverted in the subregion by explicitly specifying the region holes.
It is also possible then to refine, for example, one of the subregions.
One thing to be careful about is the application of boundary conditions like DirichletCondition with a region predicate of True. Specifying True as a predicate will apply the DirichletCondition at all boundaries, including internal boundaries. This may not necessarily be what is intended. In that case, it is better to explicitly specify the part of the boundary where the Dirichlet condition should hold. See also the Possible Issues section of DirichletCondition.
In order to create a mesh with subregions, it is sometimes useful to avoid creating region holes altogether. This can be specified by setting the "RegionHoles" option to None.
Markers
Markers are positive integer numbers that are associated with the elements in a mesh. The main purpose of markers in meshes is to detach the (boundary) predicates from actual coordinates in the PDE specification. In other words, when a PDE is specified, it can be specified in such a manner that the PDE is independent of coordinates. This functionality is useful when the simulation region is subject to change but the PDE is not.
There are three types of markers:
- Region markers are useful to specify multiple materials in a region and are part of the "MeshElements".
- Boundary markers are useful to specify NeumannValue on a region boundary and are part of the "BoundaryElements".
- Point markers are useful to specify DirichletCondition on a region boundary and are part of the "PointElements".
In this case, the triangle element has two markers, one for each element.
To illustrate the usage of material markers, a PDE with a variable diffusion coefficient is reconsidered and solved (See: Solving Partial Differential Equations with the Finite Element Method).
It is possible to specify regions markers with the "RegionMarker" option for ToElementMesh. For this, a coordinate within the region needs to be given, as well as an integer marker. Optionally, an additional maximum cell measure can be specified to refine a subregion.
Note that the lower part has a finer mesh than the upper part.
The PDE coefficient can now be specified by accessing the ElementMarker in the predicate.
Typically, specifications for boundaries and subparts in an element mesh are given by some form of a predicate. As an alternative, the boundary mesh elements can also carry markers that can be used in the boundary condition specification.
The importance of the boundary and point elements is that these are where boundary conditions are applied when solving a PDE. Generalized Neumann boundary conditions are applied by integrating over the boundary elements. Dirichlet boundary conditions are applied at the point elements.
Note that the newly inserted points and segments have the correct boundary markers set.
Note that now the PDE coefficient and boundary conditions are independent of coordinates and a new geometry can readily be put in place without the need to modify the PDE or boundary conditions.
To summarize: a set of elements is given in the form eltype[incidents,markers]. incidents is a list of the list of coordinate identities for each element. markers is a list of the same length as incidents and can be used to identify different parts of the domain or boundary where some material property might be different or a different boundary condition might apply. In evaluating on a given element, ElementMarker will effectively be replaced by the marker for that element.
Markers may be present in the "MeshElements", "BoundaryElements" and "PointElements" groups of elements. In each group, markers serve a different purpose: Markers in the "MeshElements" group of elements can be used in the PDE coefficients. Markers in the "BoundaryElements" group of elements can be used by NeumannValue or PeriodicBoundaryCondition, and markers in the "PointElements" group of elements can be used by DirichletCondition.
When specifying markers in boundary conditions, they need to be specified as equations such as 
. Several ElementMarker expressions can be combined with Boolean predicates such as And and Or. To express that a boundary condition is valid at all boundaries except a specific marker, the formulation 
can be used. Other ways of specifying ElementMarker are not supported in the current version of NDSolve.
To illustrate the process of making use of ElementMarker in a simulation, create a simple mesh that has ElementMarker in the "MeshElements", "BoundaryElements" and "PointElements".
Note that the boundary edge from node 6 to node 7 has a marker of 44, while the mesh element made up of nodes {1,6,7} also has a marker of 44 set.
As an example, a Poisson-type equation with a right-hand side that depends on the markers in the "MeshElements" is used. If a mesh element has a marker of 44, then 10 is used as the value of the right-hand side. In all other cases, 1 is used. Also a NeumannValue of 
is set on all edges that have a marker set to 44. A DirichletCondition is set on the point that has an ElementMarker of 1.
Note that the fact that both the mesh elements and boundary elements make use of a marker with number 44 does not matter. ElementMarker used in coefficients like the If statement operates on the markers present in the "MeshElements". References to ElementMarker in NeumannValue refer to markers on "BoundaryElements" and ElementMarker in DirichletCondition refers to markers in "PointElements".
Geometries may be complex. In this case, it is not always easy to identify parts of the element mesh by predicates.
By default, during boundary mesh generation an automatic grouping of contiguous boundaries is computed. The automatic grouping works by inspecting boundary normals and comparing them to each other. If the absolute value of the dot product of two normal vectors is above a threshold, which means that the vectors have close to the same direction, then the neighboring segments are considered contiguous.
Once the automatic boundary element markers are introduced, point element marks will be derived from them. For that, a node inspects the connected boundary elements' markers. If these are the same, then that marker will be used for the point element marker. At corners, a node belongs to two boundary segments with different markers. In that case, a new marker number is introduced. In previous versions, a random choice between the connected nodes and their markers was made.
This means that a DirichletCondition on the left circular part would look like this:
depending on whether the corner nodes should be considered part of the boundary condition or not.
One way to change the value of markers at the corner nodes—or any marker, for that matter—is to make use of the utility function ElementMeshResetPointElementMarker by specifying the node coordinates and the marker to be used at that coordinate.
The DirichletCondition on the left circular from above then would look like:
As a second alternative to generate a mesh where all nodes on the left edge are attributed a marker 6 and the on the right a marker 4, a "PointMakerFunction" can be used. This will change the markers based on the 
coordinate. In all other cases, the automatically generated point markers will be used.
During the generation of the full mesh, these boundary markers will then be propagated and can later be used from within NDSolve as described above.
It is possible to write functions that allow for the arbitrary placing of markers in boundary meshes. Two functions for boundary marker placement are available. One operates on boundary points, and the other operates on boundary lines in 2D and boundary faces in 3D.
Note that the order of the Which statement preceding matters. If a coordinate is on several of the boundaries, the first predicate that matches is returned.
A second function can be written to act on the boundary edges. This function gets the coordinates of the boundary element and the already-computed markers of the point function.
The union of the point markers indicates how many mesh element style directives need to be given.
Note that the marker values at the edges and points have been propagated after the full mesh was generated.
Adding mesh element markers is done by defining a "RegionMarker" option.
As an alternative, a full mesh can be generated directly.
Planar surfaces in 3D can also be subdivided and maintain their subdivision edges. For this, however, distinct markers need to be specified on the facets to indicate that the facet subdivision should be present in the full mesh.
If no markers are specified on the boundary elements, planar facets will be merged.
Note that the top surface divisions have been lost after the full mesh generation.
To maintain the surface subdivision, markers need to be added to the elements. If the markers are distinct, then the edge dividing the distinct markers will be maintained.
In this case, two of the four top surface facets have been merged because they have the same marker, 3, and surface elements with markers 1 and 2 have been propagated to the full mesh.
Sometimes it is useful to split mesh elements by their markers.
More control over the point marker computation is gained when the option "PointMarkers""BoundaryDeduced" can be specified to ToElementMesh. By specifying this method, a more rigorous approach to creating new point marker IDs is used.
Find the intersection such that only points connected to all surfaces are considered. Boundary element marker ID 7 is the large cylinder body surface.
Note that in the preceding case, the point markers on the edges that connect the cylinder and the selected surfaces are not present.
In case a constraint cannot be met, PointMarkerFromBoundaryMarker returns False, which will give a message when used in DirichletCondition: No places were found on the boundary where False was True, so DirichletCondition[u1,False] will effectively be ignored.
Element Meshes in Other Functions
Region Membership Tests
Given a point, you can test whether the point is in a region using RegionMember.
Once a region is specified, for it to be used by the finite element method, it needs to be subdivided or meshed into elements that approximately cover the region. If a mesh is already available, ElementMeshRegionMember can be used.
Test if a set of points is or is not within a meshed region.
For points that are close to the boundary of the region, wrong results are possible. This is the case when even a second-order approximation is not good enough to represent the continuous boundary properly.
One thing that can be done about this is to refine the boundary of the region.
Interpolation
When NDSolve computes a solution of a PDE via the finite element method, the returned InterpolatingFunction contains an element mesh.
It is also possible to construct an InterpolatingFunction from an ElementMesh.