Computational electromagnetics · finite element meshes

Meshing for FEM

The FEM answer is only as good as the mesh. Drag the mesh size. A real FEM solve of a microstrip line runs in your browser on each mesh. Compare the error with a converged reference.

Mesh size h / H:

Z₀ from this mesh

Error vs reference

ε_eff

Unknowns (DOF)

Mesh + 2 solves

|Z₀ error| against DOF

uniformgradedslope N^−1/2

Reference: Z₀ = , ε_eff = . Same shielded line, graded meshes up to 346,000 unknowns, extrapolated (the last step changes Z₀ by 0.003%).
Error falls slowly on a uniform mesh. The field is singular at the strip edges (like r^−1/2). Small elements there help most. Here grading reaches 2.5% error with 356 unknowns; the uniform mesh needs 4,667 for 2.9%.
P1 FEM capacitance is always too high (the discrete field energy is an upper bound). So Z₀ converges from below.

Quasi-static TEM: two Laplace solves (with and without the substrate), C and C₀ from the field energy, Z₀ = 1/(c₀√(C C₀)). Strip w = 2H, zero thickness, ε_r = 4.4, shield 12H × 8H. Linear (P1) triangles, conjugate gradients. The mesher is the Delaunay refinement of section 4 with a 20.7° angle bound.

Insert points one by one: Bowyer–Watson and Lawson flips.
March a front, then refine for quality (Ruppert).
Grade the mesh with size fields. See why slivers are a 3-D problem.
Let the solver refine, then select settings for your job.

Delaunay triangulation, one point at a time

Click to add a point. No vertex may lie inside the circumcircle of a triangle: that is the Delaunay rule (Delaunay, 1934).

Points

Triangles

Last insertion

Smallest angle

Empty-circle check (all triangles)

Bowyer (1981), Watson (1981): find the triangles whose circumcircle holds the new point. Remove them. The hole is a star-shaped cavity. Connect its boundary to the point.
Lawson (1977): split the triangle that holds the point into three. Then flip each edge that fails the circle test, until none fails.
Both give the same mesh. Of all triangulations of these points, the Delaunay one has the largest smallest angle.
Voronoi: join the circumcentres of neighbour triangles. It is the dual graph. Each cell holds the points nearest to one vertex.

Advancing front: build from the boundary in

The front starts as the boundary edges. Take the shortest front edge. Put an ideal point above it, or use a near front vertex. Add the triangle. Update the front. Repeat until the front is empty (Lo, 1985; Peraire et al., 1987).

Element size h:

Speed:

frontideal pointnew triangle

Radius ratio γ = 2r_in/R (1 = equilateral)

advancing frontDelaunay refinement (20.7°)

Front: triangles

Front: mean / min γ

Delaunay: triangles

Delaunay: mean / min γ

The front makes near-equilateral triangles in rows that follow the boundary. Its worst triangles are where two fronts meet. Delaunay refinement has a bounded worst angle but a wider spread. At this size the front wins on mean γ; at h ≥ 0.6 it does not. Gmsh “Frontal-Delaunay” (the 2-D default) mixes both: it places points like a front and connects them with Delaunay (Rebay, 1993).

Quality: Delaunay refinement with an angle bound

Delaunay triangulation alone keeps skinny triangles. Ruppert's algorithm (1995) removes them. Raise the minimum angle and watch the mesh change.

Minimum angle bound α:

0°60°

circumcircle of a bad trianglediametral circle of an encroached segment

Triangles

Smallest angle in mesh

Circumcentres added

Segment splits

Smallest angle of each triangle

Why 20.7°

ρ = R / ℓ_min = 1 / (2 sin θ_min) · ρ > √2 ⇔ θ_min < arcsin(1/(2√2)) = 20.7°

Bad triangle (smallest angle < α): insert its circumcentre. It is far from all vertices (the circle is empty), so new edges are not short.
Encroached segment (a vertex inside its diametral circle): split it at the midpoint. If a circumcentre would encroach a segment, split the segment instead.
The guarantee: split only triangles with ρ above a bound B ≥ √2. Then each new edge is longer than the edge that caused it, so the algorithm stops. That gives α up to 20.7°. Ruppert (1995) proved it for inputs with no angle under 90°; Shewchuk (2002) relaxed this to 60°. Chew's 1989 method gives 30° for uniform meshes.
In practice Triangle succeeds up to about 33° and usually does not stop above 34° (Shewchuk). Drag the bound past 34°: this demo stops at a cap of 6,000 points.

Size fields: put elements where the field changes

A size field gives the target edge length h(x, y). Gmsh takes the smallest of all size sources. Here: a Distance + Threshold field from a metal edge, and curvature sizing on a circle.

SizeMin:

SizeMax:

DistMin:

DistMax:

Curvature: elements per 2π:

Triangles

Nodes

Edge length, min / max

Mesh time

Threshold field: h against distance d to the metal edge

The same settings in the Gmsh Python API

Threshold: h = SizeMin for d ≤ DistMin, SizeMax for d ≥ DistMax, a linear ramp between (or a sigmoid with Sigmoid = 1). Gmsh computes d with the Distance field on sample points of the curve (Sampling).
Curvature: Mesh.MeshSizeFromCurvature = N asks for N elements per 2π of arc, so h = 2πR/N on the circle. The 2-D mesh then grades away from the curve.
For EM: refine at metal edges and corners (singular fields), in thin gaps, and inside high-ε_r regions where λ is short.

3-D: tetrahedra and slivers

A sliver has four good edges but almost no volume: the four vertices are nearly in one plane. Delaunay refinement in 3-D bounds the radius–edge ratio, and slivers pass that test. Drag a vertex, or flatten the tetrahedron.

Height scale (1 = as shown, 0 = flat):

Volume / regular

Dihedral min / max

Radius–edge R/ℓ_min

γ = 3r_in/R (Gmsh)

SICN (Gmsh)

κ of element K

Flatten a regular tet: metrics against height scale

γSICNR/ℓ_minκ(K) (log, right)

Why slivers hurt FEM: a dihedral angle near 0° or 180° makes the element stiffness matrix ill-conditioned. For the flattened tet, κ grows like 1/V² (its largest eigenvalue like 1/V). The gradient interpolation also gets worse: this is the maximum angle condition (Babuška & Aziz 1976; Křížek 1992). Iterative solvers then need more iterations.
No sliver problem in 2-D: a triangle with good edges cannot be flat. In 3-D, the radius–edge ratio of a flat sliver is 1/√2 = 0.71, better than TetGen's default bound of 2.0.

TetGen (Si, 2015)Delaunay refinement + optimisation

-q bounds the radius–edge ratio (default 2.0) and, as an option, the minimum dihedral angle. -O then flips edges and faces and smooths vertices; its default target is a maximum dihedral angle of 165°.

Netgen (Schöberl, 1997)advancing front + optimisation

Rule-based front placement, then edge/face swaps and smoothing. Gmsh runs it with Mesh.OptimizeNetgen = 1; Gmsh “Frontal” 3-D (4) is Netgen.

HXT (Marot et al. 2019, 2020)parallel Delaunay + optimisation

A parallel Delaunay kernel, then edge removal, Laplacian smoothing and the “growing SPR cavity” operator. Gmsh Mesh.Algorithm3D = 10.

Gmsh optimiserdefault on

Mesh.Optimize = 1 improves tetrahedra with quality below Mesh.OptimizeThreshold (default 0.3). Check with Plugin(AnalyseMeshQuality) or the SICN histogram.

Sliver exudation (Cheng et al. 2000)weighted Delaunay

Give each vertex a small weight. A weighted Delaunay mesh then has no slivers below a (very small) proven dihedral bound. A theory result; practical codes use flips and smoothing.

Adaptive refinement: let the solution choose the mesh

Solve, estimate the error of each element, mark the worst, refine, repeat. The test case is the L-shaped domain. Its re-entrant corner makes the gradient infinite, like a metal corner in EM.

Estimator

Dörfler θ:

lowhigh

Pass

DOF (free nodes)

True energy error

Estimate η / effectivity

ΔE between passes

Fitted slope (last 4)

Energy error against DOF

uniform (h/2 each pass)adaptiveN^−1/3, N^−1/2

u = r^2/3 sin(2θ/3) · |u|_H¹² = 1.8362 · P1: uniform error ∝ h^2/3 ∝ N^−1/3; adaptive ∝ N^−1/2 (optimal)

Estimate: the residual estimator sums the jumps of the normal flux across the edges of T: η_T² = ½Σ h_E²[∂_nu_h]² (Babuška & Rheinboldt 1978). ZZ averages the gradient at the nodes and measures the gap to the element gradient (Zienkiewicz & Zhu 1987). The effectivity η/error shows how honest the estimate is: near 1.0 for ZZ here, about 3.4 for the residual on the adaptive meshes (its bound has an unknown constant).
Mark (Dörfler 1996): mark the fewest elements that hold the fraction θ of the total η². Small θ: few elements per pass, a lean mesh, more passes. θ = 1 is uniform.
Refine: newest-vertex bisection (Sewell 1972; Mitchell 1991). The closure step bisects neighbours too, so there are no hanging nodes. Every element is similar to one of a few shapes of the first mesh, so the angles stay bounded. The optimal rate is proven (Binev, Dahmen & DeVore 2004; Stevenson 2007).

EM solversThe same loop in HFSS and Palace

HFSS adaptive passesstop on ΔS

Solve at one adapt frequency. Estimate the error per tetrahedron. Refine a fraction of them (Ansys training: 30% per pass by default). Stop when ΔS, the largest change of any S-parameter between two passes, is below the target (default 0.02) for the set number of converged passes. ΔS watches the output, not the field error.

Palace AMRconfig["Model"]["Refinement"]

A ZZ-type estimator on the recovered B and D fields (after Nicaise 2005). Dörfler marking with UpdateFraction (default 0.7). Stop at Tol (default 0.01) or MaxIts (default 0: no AMR). Nonconformal is true by default, with MaxNCLevels = 1; conformal refinement is for simplex meshes only.

Conformal: red–green refinement

Red: split the marked element in four. Green: bisect each neighbour to the new midpoint. Every node is a true vertex. (The demo above uses newest-vertex bisection; the result is also conforming.)

Nonconformal: hanging nodes

Only the marked element splits. The midpoints on the neighbours' edges hang (rings). The solver must constrain those DOFs to the coarse edge. It works for any element type, but the constraints add coupling.

An observed case, not a general rule

In one of Ben's Palace models, nonconformal AMR on a large mesh made the wave-port mode solve pick a spurious mode. If a port mode looks wrong after AMR, check the port field plots. On a tetrahedral mesh, also try "Nonconformal": false.

Which mesher and settings for which job?

Select a job. Settings are Gmsh option names. They are good first tries, not rules.

Job

Gmsh unstructured algorithms (Gmsh reference manual)

Option	Value	What the manual says
`Mesh.Algorithm`	1 MeshAdapt	Local edge swaps, splits and collapses. Most robust for very complex curved surfaces.
	5 Delaunay	Inserts circumcentres of the element with the largest adimensional circumradius. Fastest for very large plane surfaces; handles complex size fields well.
	6 Frontal-Delaunay	Default. Try it when element quality matters.
	7 BAMG, 8 Frontal-Delaunay for Quads, 11 Quasi-structured Quad	Anisotropic triangles; right-angle triangles for recombination into quads; quad meshes.
`Mesh.Algorithm3D`	1 Delaunay	Default. Most robust; boundary recovery by TetGen/BR; the only one for hybrid meshes with pyramids.
	4 Frontal	Netgen (Schöberl).
	10 HXT	Parallel reimplementation of Delaunay (OpenMP).
	7 MMG3D	Anisotropic tetrahedral meshes.

Embedded entities and general size fields work with 3-D Delaunay and HXT only (manual, section 1.2.1).

References

B. Delaunay, “Sur la sphère vide,” Bull. Acad. Sci. URSS, Classe Sci. Math. Nat. 6, 793–800, 1934. No DOI.
A. Bowyer, “Computing Dirichlet tessellations,” The Computer Journal 24(2), 162–166, 1981. doi:10.1093/comjnl/24.2.162
D. F. Watson, “Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes,” The Computer Journal 24(2), 167–172, 1981. doi:10.1093/comjnl/24.2.167
C. L. Lawson, “Transforming triangulations,” Discrete Mathematics 3, 365–372, 1972. doi:10.1016/0012-365X(72)90093-3; and “Software for C¹ surface interpolation,” in Mathematical Software III, Academic Press, 161–194, 1977. Max-min angle property.
S. H. Lo, “A new mesh generation scheme for arbitrary planar domains,” IJNME 21(8), 1403–1426, 1985. doi:10.1002/nme.1620210805
J. Peraire, M. Vahdati, K. Morgan, O. C. Zienkiewicz, “Adaptive remeshing for compressible flow computations,” J. Comput. Phys. 72(2), 449–466, 1987. doi:10.1016/0021-9991(87)90093-3
M. A. Yerry, M. S. Shephard, “A modified quadtree approach to finite element mesh generation,” IEEE CG&A 3(1), 39–46, 1983 doi:10.1109/MCG.1983.262997; “Automatic three-dimensional mesh generation by the modified-octree technique,” IJNME 20(11), 1965–1990, 1984 doi:10.1002/nme.1620201103. Octree meshing (not demonstrated here).
L. P. Chew, “Guaranteed-quality triangular meshes,” Tech. Rep. TR-89-983, Cornell University, 1989. No DOI.
J. Ruppert, “A Delaunay refinement algorithm for quality 2-dimensional mesh generation,” J. Algorithms 18(3), 548–585, 1995. doi:10.1006/jagm.1995.1021
J. R. Shewchuk, “Triangle: engineering a 2D quality mesh generator and Delaunay triangulator,” Applied Computational Geometry, LNCS 1148, 203–222, 1996. doi:10.1007/BFb0014497. Switch notes: cs.cmu.edu/~quake/triangle.q.html (20.7° guarantee, 33° in practice).
J. R. Shewchuk, “Delaunay refinement algorithms for triangular mesh generation,” Comput. Geom. 22(1–3), 21–74, 2002. doi:10.1016/S0925-7721(01)00047-5
J. R. Shewchuk, “Tetrahedral mesh generation by Delaunay refinement,” Proc. 14th Symp. Computational Geometry, 86–95, 1998. doi:10.1145/276884.276894 (radius–edge bound 2; slivers remain).
S.-W. Cheng, T. K. Dey, H. Edelsbrunner, M. A. Facello, S.-H. Teng, “Sliver exudation,” J. ACM 47(5), 883–904, 2000. doi:10.1145/355483.355487
H. Si, “TetGen, a Delaunay-based quality tetrahedral mesh generator,” ACM TOMS 41(2), 11:1–36, 2015. doi:10.1145/2629697. Manual: TetGen 1.5 command-line switches
J. Schöberl, “NETGEN: an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Visual. Sci. 1, 41–52, 1997. doi:10.1007/s007910050004
C. Geuzaine, J.-F. Remacle, “Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities,” IJNME 79(11), 1309–1331, 2009. doi:10.1002/nme.2579. Reference manual: gmsh.info/doc/texinfo (options and fields, checked 17 Oct 2026).
S. Rebay, “Efficient unstructured mesh generation by means of Delaunay triangulation and Bowyer–Watson algorithm,” J. Comput. Phys. 106(1), 125–138, 1993. doi:10.1006/jcph.1993.1097
C. Marot, J. Pellerin, J.-F. Remacle, “One machine, one minute, three billion tetrahedra,” IJNME 117(9), 967–990, 2019. doi:10.1002/nme.5987; C. Marot, J.-F. Remacle, “Quality tetrahedral mesh generation with HXT,” arXiv:2008.08508, 2020.
A. Johnen, C. Geuzaine, T. Toulorge, J.-F. Remacle, “Efficient computation of the minimum of shape quality measures on curvilinear finite elements,” Computer-Aided Design 103, 24–33, 2018. doi:10.1016/j.cad.2018.03.001 (SICN, SIGE). Formulas here follow the Gmsh source (qualityMeasures.cpp, qualityMeasuresJacobian.cpp).
I. Babuška, A. K. Aziz, “On the angle condition in the finite element method,” SIAM J. Numer. Anal. 13(2), 214–226, 1976. doi:10.1137/0713021; M. Křížek, “On the maximum angle condition for linear tetrahedral elements,” SIAM J. Numer. Anal. 29(2), 513–520, 1992. doi:10.1137/0729031
I. Babuška, W. C. Rheinboldt, “Error estimates for adaptive finite element computations,” SIAM J. Numer. Anal. 15(4), 736–754, 1978. doi:10.1137/0715049
O. C. Zienkiewicz, J. Z. Zhu, “A simple error estimator and adaptive procedure for practical engineering analysis,” IJNME 24(2), 337–357, 1987. doi:10.1002/nme.1620240206; “The superconvergent patch recovery and a posteriori error estimates, Part 1,” IJNME 33(7), 1331–1364, 1992. doi:10.1002/nme.1620330702
W. Dörfler, “A convergent adaptive algorithm for Poisson's equation,” SIAM J. Numer. Anal. 33(3), 1106–1124, 1996. doi:10.1137/0733054
W. F. Mitchell, “Adaptive refinement for arbitrary finite-element spaces with hierarchical bases,” J. Comput. Appl. Math. 36(1), 65–78, 1991. doi:10.1016/0377-0427(91)90226-A. E. G. Sewell, PhD thesis, Purdue University, 1972 (newest-vertex bisection; not checked first-hand).
M.-C. Rivara, “Algorithms for refining triangular grids suitable for adaptive and multigrid techniques,” IJNME 20(4), 745–756, 1984. doi:10.1002/nme.1620200412 (longest-edge bisection).
P. Binev, W. Dahmen, R. DeVore, “Adaptive finite element methods with convergence rates,” Numer. Math. 97, 219–268, 2004. doi:10.1007/s00211-003-0492-7; R. Stevenson, “Optimality of a standard adaptive finite element method,” Found. Comput. Math. 7, 245–269, 2007. doi:10.1007/s10208-005-0183-0
W. F. Mitchell, “A collection of 2D elliptic problems for testing adaptive grid refinement algorithms,” Appl. Math. Comput. 220, 350–364, 2013. doi:10.1016/j.amc.2013.05.068 (the L-shaped benchmark).
S. Nicaise, “On Zienkiewicz–Zhu error estimators for Maxwell's equations,” C. R. Math. 340, 697–702, 2005. doi:10.1016/j.crma.2005.03.016
Palace documentation: awslabs.github.io/palace (reference: error estimation and AMR; config["Model"]["Refinement"] defaults from the configuration schema, checked 17 Oct 2026).
Ansys HFSS Getting Started, lecture 3 “FEM and mesh solution” (2020 R2), Ansys Innovation Space. ΔS default 0.02 and 30% refinement per pass; vendor training slides, details vary by release.
F. Hecht, BAMG: bidimensional anisotropic mesh generator, INRIA, 1998; C. Dobrzynski, MMG3D user guide, INRIA, 2012. As cited by the Gmsh manual.