Element Mesh Generation

In order to use mesh generation functionality, the finite element method (FEM) package needs to be loaded.

Load the package.

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:

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.

Set up a region.
Set up a PDE operator.
Specify boundary conditions.
Solve the PDE.
Plot a contour plot of the solution with the element mesh from the interpolation function on top.

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.

Extract the ElementMesh from an interpolating function.

Instead of specifying an implicit parametric region, it is also possible to specify an explicit ElementMesh. This can be done by using ToElementMesh.

Define an explicit mesh with four triangle elements.
Show a wireframe of the element mesh.

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.

Solve the PDE with an ElementMesh.
This makes a contour plot of the solution and plots the element mesh.

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.

Solve the PDE with options given to influence the element mesh generation.

As an alternative, the element mesh can be generated prior to the simulation and given to NDSolve.

Generate an element mesh, visualize the mesh, and solve the PDE.
Solve the time-dependent PDE with options given to influence the element mesh generation.

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.

Investigate the mesh order of an ElementMesh.

The ability to deal with higher-order elements, expressed through the mesh order, has two distinct advantages:

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:

Line meshes

For a 1D mesh, the mesh elements ei are LineElement. The boundary elements bi are then PointElement.

Create a line element mesh with three coordinates and two line elements.
Visualize the element mesh wireframe.

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.

Create a triangle element mesh with six triangle elements.
Inspect the mesh order.
Visualize the element mesh wireframe and the coordinates with their incident numbers.

The element-element connectivity provides neighboring information about elements.

Visualize the element mesh wireframe with the element's IDs in black and node IDs in red.

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.

Extract the element connectivity information.

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.

Create a triangle element mesh with six triangle elements.
Visualize the element mesh wireframe with the element's markers in red.

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.

Specify coordinates for a second-order mesh.
Specify six quadratic triangle elements without markers.
Create a second-order mesh.
Visualize the second-order mesh.

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.

Construct coordinates on a circular perimeter.
The corresponding incidents.
Create a quad element mesh.
Visualize the element mesh wireframe.

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.

Create a mixed-element mesh.
Visualize the element mesh wireframe.

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.

Extract the element incidents from the previous elements and add integer 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.

Extract the element connectivity.

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.

A boundary mesh with six coordinates and a bounding line.
Visualize the boundary mesh.
Convert the boundary mesh to a full mesh and visualize.
A boundary mesh with six coordinates, a bounding line and a separation.
Visualize the boundary mesh.
Convert the boundary mesh to a full mesh, respecting the separation and visualizing.

For a boundary mesh, you can add arbitrary points inside the bounding region to be part of the mesh.

Create boundary coordinates and some points placed on a circle.
Create the boundary mesh.
All coordinates given are referenced point elements.
Visualize the boundary mesh and the point elements.
Visualize the full mesh.

Tetrahedron meshes

Three-dimensional manual mesh creation follows the same ideas as in one and two dimensions.

A linear tetrahedron element mesh with markers.
Visualizing the index of the coordinates at their respective positions.
Create the mesh.
Visualize the mesh with the elements' markers.

Hexahedron meshes

Construct coordinates on a circular perimeter.
The linear incidents.
Create a hexahedron element mesh.
Visualize the boundary element mesh wireframe.

Mixed-element type meshes in 3D

For a 3D mesh, the mesh elements ei can be either a TetrahedronElement or HexahedronElement but not both unless a PrismElement is used to connect them.

Boundary meshes in 3D

Boundary meshes are useful to generate full meshes. In 3D, the boundary elements bi are TriangleElement or QuadElement.

Here are coordinates in 3D space.
Quad elements as boundary elements.
Create the boundary mesh.
Visualize the full mesh.

Symbolic Regions

The usage of symbolic region representations and their conversion to an ElementMesh are straightforward. Any graphics primitive can be used to express a symbolic region.

Create a symbolic region of a Disk.
Convert the symbolic region to an ElementMesh.
Visualize the mesh.
Create and visualize a mesh from a 3D symbolic region.

Boolean combinations of symbolic regions are also possible.

Create and visualize a mesh from a region difference.

Please consult the reference pages of ToBoundaryMesh and ToElementMesh for many more examples.

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.

Specify a symbolic region.
Generate a mesh.
Inspect the full mesh.

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.

Use a higher PrecisionGoal to resolve the small feature.
Inspect the full mesh.

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.

Generate a NumericalRegion from the symbolic region representation.

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.

Generate boundary meshes of the separate regions.
Extract the boundary elements from the boundary meshes.

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.

Manually generate a boundary mesh.
Visualize the boundary mesh and a zoomed region of the boundary mesh.
Engage the numerical region with the new boundary mesh.
Inspect the symbolic region in the numerical region.
Inspect the boundary element mesh in the numerical region.
Generate the full mesh from the numerical region. Even though the boundary elements were first order, the elements can be made second order because the numerical region allows correct placement of the requisite mid-side nodes.

Note that there are far fewer mesh elements for the same approximation quality.

Inspect the wireframe mesh.
Verify that the mesh is second-order mesh.

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.

Special Purpose Meshes

Region products

Another alternative to create 2D and 3D meshes is by making use of region products. If a 1D mesh is multiplied by another 1D mesh a 2D mesh results and if a 2D mesh is multiplied by a 1D mesh a 3D mesh results.

Create 1D mesh.
Form the region product of two 1D meshes.
Visualize the 2D region product.

Region products can also be made from symbolic regions.

Create and visualize a coarse mesh of an annulus.
Create the region product of a 2D mesh with a coarse 1D mesh.

Graded meshes

The function ToGradedMesh generates 1D graded meshes. A graded mesh is a mesh that has a non uniform distribution of nodes.

Create a graded mesh with a high coordinate concentration on the left.
Visualize the mesh.

Note the higher concentration of nodes on the left. Other forms of alignment, like "Right", "BothEnds", "Central" and "Uniform" are possible. Also custom alignment functions can be give.

Graded meshes can be useful in a number of scenarios. To mimic an infinite domain one can model a relatively large region and only have a high concentration of nodes in one area. This is shown in an example on the reference page of ToGradedMesh. Another area where graded meshes are useful is if a material boundary has a discontinuity. When a graded mesh with a high concentration of nodes at the problematic region is generated then the approximation quality of the solution can be improved, sometimes substantially.

If is also possible to create 2D and 3D graded meshes with the help of ElementMeshRegionProduct.

Create a 1D graded mesh with a "Central" point distribution.
Create a second 1D graded mesh with a "Central" point distribution and 50 elements.
Construct a region product.
Visualize the region product of the anisotropic mesh.
Create a third 1D graded mesh with a "BothEnds" point distribution, 50 elements and an endpoint distance of 1/200.
Create the region product.
Visualize the 3D graded mesh.

For more application examples please see the applications section of ToGradedMesh and ElementMeshRegionProduct.

Meshes for perfectly matched layers

A perfectly matched layer (PML) extends a mesh and a corresponding partial differential equation to allow for an attenuation of a signal in the perfectly matched layer. This is useful to mimic domains of infinite extend but still use a finite sized domain. The creation of these type of meshes are best explained by looking at specific examples. For acoustics there are explanations of the usage and creation of PMLs in the time domain and PMLs in the frequency domain. There is also a 3D PML electromagnetics example.

Meshes from contours

Sometimes a mesh is needed that has the contour or an object as its boundary. Such meshes are typically used in fluid simulations, for example, to simulate the flow around an object, like a plane. This section demonstrates how a mesh around the island of Pico can be created. The data is retrieved from the geo objects of the island.

Specify the position of the island and visualize the geo graphics.
Extract the geo elevation data from the geo graphics.
Clip and scale the data.
Downsample the data.
Visualize the data with a "ReliefMap" texture.
Create a discretization of the relief plot.
Create the sides and top for a bounding box.

Note that the box is open at the bottom.

Create the union of the bounding box and the relief discretization.
Create the mesh.
Visualize the final mesh.

Next, a visualization is generated that shows the relief plot with a part of the mesh removed. For this, you compute the center of mass of each element and select only those elements that you would like to see in the final visualization.

Compute the center of mass of each element.
Choose which elements should be visible in the final visualization.
Show the partial visualization of the mesh and the relief plot.

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.

Numerically integrate over a Rectangle.
Sum the area of an ElementMesh discretization of a Rectangle.

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 difference of the exact value and the sum of the area of an ElementMesh discretization of a Disk.

The discrepancy is there, despite that the boundary points are close on the boundary.

Inspect the approximation quality of the boundary coordinates.
Inspect the influence the mesh order has on the approximation of a Disk.

In the case of a Disk, the mesh order plays a crucial role, since the ElementMesh can hold curved boundary mesh elements.

Inspect the influence the mesh granularity on the boundary has on the approximation of a Disk.

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.

Specify an AccuracyGoal for a first-order approximation of a Disk.
Inspect the influence of an AccuracyGoal for a second-order approximation of a Disk.

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.

Inspect the influence of "ImproveBoundaryPositions" for a second-order approximation of a Disk.

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.

Set up a region with cusps.
Use the "Continuation" boundary mesh generator and visualize the mesh.
Compare the area of the region Ω with the approximation.

As a contrast, compare the result to the fast region plot method.

Use the "RegionPlot" boundary mesh generator and visualize the mesh.

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).

Use the default boundary mesh generator with more sample points and visualize the mesh.
Compare the area of the region Ω with the approximation.

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]:

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.

The optimal triangle element.
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.

The optimal quad element.
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.

The optimal tetrahedron.
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.

The optimal hexahedron element.
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. 309319.

[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.

Convert a space shuttle's orbiter to an element mesh.
Visualize the element mesh surface.
Compute the element mesh quality and overview quantities such as the minimum quality.
Visualize the mesh element quality distribution.

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.

Improve the element mesh quality by requesting a smaller minimum edge ratio. This will result in more elements.
Compute quality data and overview quality quantities.

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.

Visualize the boundary mesh.

To visualize the elements that are below a certain threshold, the poor-quality elements are selected from the mesh element quality list and visualized.

Visualize the mesh elements that have a quality of less than 0.1.
Visualize the mesh elements that have a quality of less than 0.1 combined with all mesh elements.

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 intersect the internal material boundaries. To illustrate this, a PDE with a variable diffusion coefficient is reconsidered and solved over two material regions. (See: Solving Partial Differential Equations with the Finite Element Method.) One region respects the internal boundary, while the other does not.

Create and visualize a boundary element mesh with an internal boundary and without an internal boundary.

The diffusion coefficient has a jump discontinuity at .

Set up a diffusion coefficient that is space dependent.
Create and visualize the element meshes.
Specify a PDE operator and boundary conditions.
Solve the equation over each mesh.
Visualize the solution.
Visualize the difference between the two solutions.

The next example shows a circular region with a subregion and holes inside the subregion.

Set up the outer annulus and the inner region holes.
Display the region.

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.

Specify the region as an implicit region and create an element mesh.

In the next step, what is a region hole and what is not is inverted in the subregion by explicitly specifying the region holes.

Create a mesh with specified region holes.

It is also possible then to refine, for example, one of the subregions.

Refine a subregion.

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.

Avoid creating region holes and visualize the mesh element markers.

Additional information about how to create material regions in three dimensions can be found in the OpenCascadeLink.

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:

Evaluating Expressions with Markers

Sometimes one would like to evaluate expressions that depend on ElementMarker. This can be done with the function EvaluateOnElementMesh.

Create a mesh with markers.
Visualize the mesh with markers.
Create an expression with ElementMarker.
Evaluate the expression over the mesh.

Evaluating the expression on the mesh returns an interpolating function that can be used like any other function.

Numerically integrate over the function.
Check that the integration is correct.

The returned interpolating function behaves the same way as other interpolating functions.

Evaluate the interpolating function outside of the data range.

This behavior can be changed as described in ElementMeshInterpolation.

Evaluate the expression over the mesh with a specified extrapolation handler.
Evaluate the interpolating function outside of the data range.

Note that the usage of EvaluateOnElementMesh can lead to overshooting and undershooting, as explained in the section Overshoot/Undershoot Issue for Discontinuous Coefficients.

Units

Sometimes a mesh needs to be rescaled. For example, an imported mesh is not in the correct scale or, to avoid tolerance issues, a mesh is generated in a different scale and needs to be adjusted in a subsequent step.

Generate a boundary mesh.
Inspect the length unit of the mesh.
Set a length unit.
Rescale the coordinates of the mesh.

Keep in mind, though, that for a subsequent finite element analysis it might be advisable to rescale the PDE coefficients for better numerical stability of the results.

Element Meshes in Other Functions

Region Membership Tests

Given a point, you can test whether the point is in a region using RegionMember.

Create a region and test if the point is in the region .

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.

Create a mesh and test if the point is in the mesh.

Test if a set of points is or is not within a meshed region.

Test if a number of points are inside 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.

Test a large number of points for region inclusion and extract those that differ from the inclusion of the exact region.
Visualize the mesh with the points close to the boundary.

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.

Solve a PDE.
This gets the ElementMesh from the InterpolatingFunction returned for the Poisson equation.

It is also possible to construct an InterpolatingFunction from an ElementMesh.

Create an ElementMesh.
Evaluate a function over the coordinates of the mesh.
Visualize the InterpolatingFunction.

Interpolating functions can be used as coefficients in PDEs. When the mesh in the interpolating function is the same as the mesh for the PDE solution process then a solution process is faster.

Measure the timing used to solve a PDE where the same mesh is used for the PDE and the interpolating function.
Measure the timing used to solve a PDE where the same mesh is not used for the PDE and the interpolating function.

You can see that using the same mesh for an interpolating function and a PDE solution is more efficient.