Element Mesh Generation

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

Load the package.
In[1]:=
Click for copyable input

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 way 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 mesh creation
  • explicit creation

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.

Set up a region.
In[3]:=
Click for copyable input
Out[3]=
Set up a PDE operator.
In[4]:=
Click for copyable input
Out[4]=
Specify boundary conditions.
In[5]:=
Click for copyable input
Out[5]=
Solve the PDE.
In[6]:=
Click for copyable input
Out[6]=
Plot a contour plot of the solution with the element mesh from the interpolation function on top.
In[7]:=
Click for copyable input
Out[7]=

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.
In[8]:=
Click for copyable input
Out[8]=

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.
In[9]:=
Click for copyable input
Out[9]=
Show a wireframe of the element mesh.
In[10]:=
Click for copyable input
Out[10]=

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 the explicit mesh only defined.

Solve the PDE with an ElementMesh.
In[11]:=
Click for copyable input
This makes a contour plot of the solution and plots the element mesh.
In[12]:=
Click for copyable input
Out[12]=

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

All options for ToElementMesh and ToBoundaryMesh can be given to NDSolve directly.

Solve the PDE with options given to influence the element mesh generation.
In[13]:=
Click for copyable input

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.
In[14]:=
Click for copyable input
Out[15]=
In[16]:=
Click for copyable input
Out[16]=

Comparing ElementMesh and MeshRegion

Before a deeper discussion of the ElementMesh functionality is given, it is instructional to compare the ElementMesh object to mesh regions.

The default for an ElementMesh is 2.

Investigate the mesh order of an ElementMesh.
In[17]:=
Click for copyable input
Out[17]=

The ability to deal with higher-order elements has two distinct advantages:

  • a higher-quality numerical approximation
  • the ability to better approximate curved boundaries
Compare an analytical solution of a PDE with a numerical solution computed with a first-order mesh and a second-order mesh.
In[18]:=
Click for copyable input
Out[22]=
Compute the area of an approximation of a Disk based on a MeshRegion and an ElementMesh.
In[23]:=
Click for copyable input
Out[23]=
In[24]:=
Click for copyable input
Out[24]=

The ElementMesh is able to approximate the Disk better because it uses a second-order approximation.

The conversion between MeshRegion and ElementMesh is easy.

Convert a MeshRegion to an ElementMesh.
In[25]:=
Click for copyable input
Out[25]=
Convert an ElementMesh to a MeshRegion.
In[26]:=
Click for copyable input
Out[26]=

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.

Create a boundary and a full ElementMesh.
In[31]:=
Click for copyable input
Out[31]=
In[32]:=
Click for copyable input
Out[32]=
Inspect the mesh elements of a boundary and a full mesh.
In[33]:=
Click for copyable input
Out[33]=
In[34]:=
Click for copyable input
Out[34]//Short=

A boundary ElementMesh does not need to be a closed curve.

Create a boundary ElementMesh that is not a closed curve.
In[35]:=
Click for copyable input
Out[35]=

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.

Create a MeshRegion with a lower-dimensional component.
In[36]:=
Click for copyable input
Out[36]=
Inspect the dimensional components for a MeshRegion.
In[37]:=
Click for copyable input
Out[37]=
Convert and visualize a MeshRegion with lower-dimensional components to an ElementMesh.
In[38]:=
Click for copyable input
Out[38]=

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.

Create and visualize a boundary ElementMesh with internal structure.
In[39]:=
Click for copyable input
Out[40]=
Create and visualize a full ElementMesh with an internal structure.
In[41]:=
Click for copyable input
Out[42]=

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 use an ElementMesh. Some of the remaining items will be further discussed in this tutorial.

Pros

  • memory-efficient storage
  • support for higher-order mesh elements
  • support for curved mesh elements
  • if needed, a conversion to e.g MeshRegion is easily possible
  • works for both boundary and full meshes
  • setting markers is easy
  • influence the node ordering

Cons

  • not in system context; a package needs to be loaded
  • only usable within the finite element method framework

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:

An element is specified by its type, for example TriangleElement and a list of integer lists.

In[43]:=
Click for copyable input
Out[43]=

The integers are indices that correspond to coordinates. In the above case, the TriangleElement contains two triangle elements, one made up of the incidents , and one made up of the incidents .

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[44]:=
Click for copyable input
Out[44]=

The element indices (also called incidents) correspond to the coordinates. The incidents of the second triangle element are and refer to coordinates , , and , respectively.

Create an element mesh.
In[45]:=
Click for copyable input
Out[45]=

The concept of the element incidents is closely related to that of a GraphicsComplex.

Visualize the element mesh wireframe.
In[46]:=
Click for copyable input
Out[46]=

Line Meshes

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

Create a line element mesh with three coordinates and two line elements.
In[47]:=
Click for copyable input
Out[47]=
Visualize the element mesh wireframe.
In[48]:=
Click for copyable input
Out[48]=

Triangle Meshes

For a 2D mesh, the mesh elements are TriangleElement or QuadElement. The boundary elements 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.
In[49]:=
Click for copyable input
Out[49]=
Inspect the mesh order.
In[50]:=
Click for copyable input
Out[50]=
Visualize the element mesh wireframe and the coordinates with their incident numbers.
In[51]:=
Click for copyable input
Out[51]=

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.
In[52]:=
Click for copyable input
Out[52]=
Visualize the element mesh wireframe with the element's markers in red.
In[53]:=
Click for copyable input
Out[53]=

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.
In[54]:=
Click for copyable input
Specify six quadratic triangle elements without markers.
In[55]:=
Click for copyable input
Create a second-order mesh.
In[56]:=
Click for copyable input
Out[56]=
Visualize the second-order mesh.
In[57]:=
Click for copyable input
Out[57]=

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.
In[58]:=
Click for copyable input
The corresponding incidents.
In[60]:=
Click for copyable input
Create a triangle element mesh.
In[61]:=
Click for copyable input
Out[61]=
Visualize the element mesh wireframe.
In[62]:=
Click for copyable input
Out[62]=

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 can be a combination of TriangleElement and QuadElement. The boundary elements are then LineElement.

Create a mixed element mesh.
In[63]:=
Click for copyable input
In[64]:=
Click for copyable input
In[66]:=
Click for copyable input
Out[66]=
Visualize the element mesh wireframe.
In[67]:=
Click for copyable input
Out[67]=

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.
In[68]:=
Click for copyable input
In[69]:=
Click for copyable input
Out[69]=
In[70]:=
Click for copyable input
Out[70]=

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[71]:=
Click for copyable input
Out[71]=

In the above 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.
In[72]:=
Click for copyable input
Out[72]=

Boundary Meshes in 2D

Boundary meshes are useful to generate full meshes. In 2D, the boundary elements are the LineElement.

A boundary mesh with six coordinates and a bounding line.
In[73]:=
Click for copyable input
Out[73]=
Visualize the boundary mesh.
In[74]:=
Click for copyable input
Out[74]=
Convert the boundary mesh to a full mesh and visualize.
In[75]:=
Click for copyable input
Out[75]=
A boundary mesh with six coordinates, a bounding line, and a separation.
In[76]:=
Click for copyable input
Out[76]=
Visualize the boundary mesh.
In[77]:=
Click for copyable input
Out[77]=
Convert the boundary mesh to a full mesh, respecting the separation and visualizing.
In[78]:=
Click for copyable input
Out[78]=

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.
In[79]:=
Click for copyable input
Create the boundary mesh.
In[80]:=
Click for copyable input
Out[80]=
All coordinates given are referenced point elements.
In[81]:=
Click for copyable input
Out[81]=
Visualize the boundary mesh and the point elements.
In[82]:=
Click for copyable input
Out[82]=
Visualize the full mesh.
In[83]:=
Click for copyable input
Out[83]=

Tetrahedron Meshes

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

A linear tetrahedron element mesh with markers.
In[84]:=
Click for copyable input
Visualizing the index of the coordinates at their respective positions.
In[85]:=
Click for copyable input
Out[85]=
Create the mesh.
In[86]:=
Click for copyable input
Out[86]=
Visualize the mesh with the elements' markers.
In[87]:=
Click for copyable input
Out[87]=

Hexahedron Meshes

Construct coordinates on a circular perimeter.
In[88]:=
Click for copyable input
The linear incidents.
In[90]:=
Click for copyable input
Create a hexahedron element mesh.
In[91]:=
Click for copyable input
Out[91]=
Visualize the boundary element mesh wireframe.
In[92]:=
Click for copyable input
Out[92]=

Mixed Element Type Meshes in 3D

For a 3D mesh, the mesh elements 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 are TriangleElement or QuadElement.

Here are coordinates in 3D space.
In[93]:=
Click for copyable input
Out[93]=
Quad elements as boundary elements.
In[94]:=
Click for copyable input
Out[94]=
Create the boundary mesh.
In[95]:=
Click for copyable input
Out[95]=
In[96]:=
Click for copyable input
Out[96]=
Visualize the full mesh.
In[97]:=
Click for copyable input
Out[97]=

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.
In[98]:=
Click for copyable input
Out[98]=
Sum the area of an ElementMesh discretization of a Rectangle.
In[99]:=
Click for copyable input
Out[99]=

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.
In[100]:=
Click for copyable input
Out[100]=

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

Inspect the approximation quality of the boundary coordinates.
In[101]:=
Click for copyable input
Out[102]=
Inspect the influence the mesh order has on the approximation of a Disk.
In[103]:=
Click for copyable input
Out[103]=

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.
In[104]:=
Click for copyable input
Out[104]=

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.
In[105]:=
Click for copyable input
Out[105]=
Inspect the influence of an AccuracyGoal for a second-order approximation of a Disk.
In[106]:=
Click for copyable input
Out[106]=

Internally, ToBoundaryMesh is called first, 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 for a second-order approximation of a Disk.
In[107]:=
Click for copyable input
Out[107]=

In this case, the result is not different from a first-order approximation.

For the creation of the boundary mesh, two boundary mesh generators are available. The default is based on RegionPlot and called . This provides a fast boundary approximation. The second boundary mesh generator is called 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.
In[108]:=
Click for copyable input
Out[108]=
Use the boundary mesh generator and visualize the mesh.
In[109]:=
Click for copyable input
Out[110]=
Compare the area of the region Ω with the approximation.
In[111]:=
Click for copyable input
Out[111]=

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

Use the boundary mesh generator and visualize the mesh.
In[112]:=
Click for copyable input
Out[113]=

The visual inspection already shows that the cusps are not well resolved. The 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.
In[114]:=
Click for copyable input
Out[115]=
Compare the area of the region Ω with the approximation.
In[116]:=
Click for copyable input
Out[116]=

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

  • large angles (near 180°) can cause large interpolation errors
  • small angles (near 0°) can cause the finite element stiffness matrix to be ill conditioned
  • smaller elements offer more accuracy but cost more computational time
  • small or skinny elements can induce instabilities in explicit time integration methods

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.

Construct a mesh and visualize a mesh.
In[117]:=
Click for copyable input
In[120]:=
Click for copyable input
Out[120]=
In[121]:=
Click for copyable input
Out[121]=
Compute the quality of each element in an element mesh.
In[122]:=
Click for copyable input
Out[122]//Short=

The result of a computation is a quality estimate for each mesh element grouped by mesh elements.

The mean and standard deviation of the element mesh quality can be computed.
In[123]:=
Click for copyable input
Out[123]=
Compute the minimum and maximum element mesh quality.
In[124]:=
Click for copyable input
Out[124]=
Visualize the mesh element quality in a histogram.
In[125]:=
Click for copyable input
Out[125]=

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.

The optimal triangle element.
In[126]:=
Click for copyable input
Out[127]=
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.
In[128]:=
Click for copyable input
Out[129]=
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.
In[130]:=
Click for copyable input
Out[131]=
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.
In[132]:=
Click for copyable input
Out[133]=
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.
In[134]:=
Click for copyable input
Out[135]=
Visualize the element mesh surface.
In[136]:=
Click for copyable input
Out[136]=
Compute the element mesh quality and overview quantities such as the minimum quality.
In[137]:=
Click for copyable input
Out[138]=
Visualize the mesh element quality distribution.
In[139]:=
Click for copyable input
Out[139]=

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.
In[140]:=
Click for copyable input
Out[140]=
Compute quality data and overview quality quantities.
In[141]:=
Click for copyable input
Out[142]=

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.

In[143]:=
Click for copyable input
Out[143]=
Visualize the boundary mesh.
In[144]:=
Click for copyable input
Out[144]=

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.
In[145]:=
Click for copyable input
Out[146]=
Visualize the mesh elements that have a quality of less than 0.1 combined with all mesh elements.
In[147]:=
Click for copyable input
Out[148]=

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.

Create and visualize a boundary element mesh with an internal boundary and without an internal boundary.
In[149]:=
Click for copyable input
Out[154]=

The diffusion coefficient has a jump discontinuity at .

Set up a diffusion coefficient that is space dependent.
In[155]:=
Click for copyable input
Out[155]=
Create and visualize the element meshes.
In[156]:=
Click for copyable input
Out[158]=
Specify a PDE operator and boundary conditions.
In[159]:=
Click for copyable input
Solve the equation over each mesh.
In[161]:=
Click for copyable input
In[162]:=
Click for copyable input
Visualize the solution.
In[163]:=
Click for copyable input
Out[163]=
Visualize the difference between the two solutions.
In[164]:=
Click for copyable input
Out[164]=

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 .
In[165]:=
Click for copyable input
Out[166]=

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, can be used.

Create a mesh and test if the point is in the mesh.
In[167]:=
Click for copyable input
Out[168]=

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.
In[169]:=
Click for copyable input
Out[173]=

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.
In[174]:=
Click for copyable input
Out[177]=
Out[179]=
Visualize the mesh with the points close to the boundary.
In[180]:=
Click for copyable input
Out[180]=

One thing that can be done about this is to refine the boundary of the region.

In[181]:=
Click for copyable input
Out[181]=
In[182]:=
Click for copyable input
Out[182]=
In[183]:=
Click for copyable input
Out[184]=
Out[186]=

Interpolation

When NDSolve computes a solution of a PDE via the finite element method, the returned InterpolatingFunction contains an element mesh.

Solve a PDE.
In[187]:=
Click for copyable input
This gets the ElementMesh from the InterpolatingFunction returned for the Poisson equation.
In[188]:=
Click for copyable input
Out[188]=

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

Create an ElementMesh.
In[189]:=
Click for copyable input
Out[189]=
Evaluate a function over the coordinates of the mesh.
In[190]:=
Click for copyable input
In[191]:=
Click for copyable input
Out[191]=
Visualize the InterpolatingFunction.
In[192]:=
Click for copyable input
Out[192]=

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.

Create a triangle mesh with two element and markers.
In[193]:=
Click for copyable input
Out[193]=

In this case, the triangle element has two markers, one for each element.

Visualize the mesh and the markers.
In[194]:=
Click for copyable input
Out[194]=

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

Create and visualize a boundary element mesh with an internal boundary.
In[195]:=
Click for copyable input
Out[197]=

It is possible to specify regions markers with the 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.

Create and visualize an element mesh with internal markers.
In[198]:=
Click for copyable input
Out[199]=

Note that the lower part has a finer mesh than the upper part.

Extract the markers from the mesh elements with .
In[200]:=
Click for copyable input
Out[200]=

The PDE coefficient can now be specified by accessing the in the predicate.

Specify the PDE coefficient with an .
In[201]:=
Click for copyable input
Out[201]=
Specify a PDE operator and boundary conditions.
In[202]:=
Click for copyable input
Solve the equation.
In[204]:=
Click for copyable input
Visualize the solution.
In[205]:=
Click for copyable input
Out[205]=

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.

Create and visualize a boundary element mesh with an internal boundary and boundary markers.
In[206]:=
Click for copyable input
Out[208]=
Create and visualize an element mesh with internal markers.
In[209]:=
Click for copyable input
Out[210]=
Visualize the boundary and its markers of the mesh.
In[211]:=
Click for copyable input
Out[211]=
Extract the markers from the boundary elements with .
In[212]:=
Click for copyable input
Out[212]=
Extract the markers from the point elements with .
In[213]:=
Click for copyable input
Out[213]=

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.

Specify the PDE coefficient with an ElementMarker.
In[214]:=
Click for copyable input
Out[214]=
Specify a PDE operator.
In[215]:=
Click for copyable input
Out[215]=
Specify boundary conditions with markers.
In[216]:=
Click for copyable input
Solve the equation.
In[217]:=
Click for copyable input
Out[217]=
Visualize the solution.
In[218]:=
Click for copyable input
Out[218]=

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

Geometries may be complex. In this case, it is not always easy to identify parts of the element mesh by predicates.

Create a list of boundaries and visualize them.
In[219]:=
Click for copyable input
Out[220]=
Create a boundary mesh by setting the boundaries smaller than 0 and combine them with And.
In[221]:=
Click for copyable input
Out[221]=
Visualize the boundary mesh.
In[222]:=
Click for copyable input
Out[222]=

It is possible to write functions that allow for placing of markers in boundary meshes. During the generation of the full mesh, these boundary markers will then be propagates and can later be used from within NDSolve as described above. 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.

Write a function that, given coordinates of the boundary mesh, will return integer markers if the coordinates satisfy a predicate.
In[223]:=
Click for copyable input
Apply the point marker function to the coordinates of the boundary mesh.
In[224]:=
Click for copyable input
Out[224]=
Visualize the integer markers attributed to the coordinates of the boundary mesh.
In[225]:=
Click for copyable input
Out[225]=

Note that the order of the Which statement above 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.

Here is a function that computes integer markers for edges depending only on the integer markers of the boundary points.
In[226]:=
Click for copyable input
The point and edge marker functions can be given to ToBoundaryMesh.
In[227]:=
Click for copyable input
Out[227]=
Inspecting the boundary elements of the mesh reveals the boundary markers.
In[228]:=
Click for copyable input
Out[228]=
Visualize the boundary mesh with its markers.
In[229]:=
Click for copyable input
Out[229]=
Inspecting the point elements of the mesh reveals the point markers.
In[230]:=
Click for copyable input
Out[230]=
Visualize the boundary mesh with its markers.
In[231]:=
Click for copyable input
Out[231]=
Convert the boundary mesh to a full mesh.
In[232]:=
Click for copyable input
Out[232]=
Inspect the union of the point markers.
In[233]:=
Click for copyable input
Out[233]=
Visualize the element mesh and the boundary markers.
In[234]:=
Click for copyable input
Out[234]=
Inspect the union of the point markers.
In[235]:=
Click for copyable input
Out[235]=

The union of the point markers indicates how many mesh element style directives need to be given.

Visualize the element mesh and the point markers.
In[236]:=
Click for copyable input
Out[236]=

Note that the marker values at the edges and points have been propagated after the full mesh was generated.

As an alternative, a full mesh can be generated directly.

Directly create a full mesh with point and boundary markers.
In[237]:=
Click for copyable input
Out[237]=