Finite Element Method Magnetics

A 2D constrained Delaunay triangulation and quality mesh generation tool, implemented in C++17. File format compatible with Shewchuk's Triangle with extensions for arc segments, local feature size constraints, and periodic boundary conditions. Also reads FEMM problem files (.fem, .fee, .feh, .fec) directly.

A Python implementation (tangle.py) is also included with identical functionality (.poly/.node input only).

• David Meeker — author
• Claude Code — development assistance

Tangle is backward-compatible with Triangle's .node, .poly, .ele, .edge, and .neigh file formats, for compatibility with existing workflows. Tangle also provides additional capabilities that are practically important to formulation of meshes for the solution of finite element problems (particularly in the area of electric machines):

• Local feature size (LFS) constraints — per-vertex and per-segment mesh density control, so that refinement can be concentrated where it matters without over-meshing the rest of the domain.
• Periodic boundary conditions (PBCs) — paired boundary chains with synchronized segment splitting during refinement, producing matching node-to-node correspondence for periodic or anti-periodic coupling.
• Air gap elements (AGEs) — structured element layers across motor air gaps with uniformly-spaced chord segments that are protected from refinement splitting.
• Arc segments — circular arcs defined by two endpoints and a subtended angle, automatically discretized into chord segments.
• Native FEMM file reading — directly meshes .fem, .fee, .feh, and .fec files, handling arc discretization, SmartMesh, boundary properties, and block labels without a separate preprocessor.

These extensions are described in polyformat.txt, which documents the full extended .poly format.

g++ -O3 -std=c++17 -static -o tangle tangle.cpp float256.cpp -lm

No external dependencies. Requires a C++17 compiler (GCC, Clang, or MSVC).

Source code repository at github, pre-built Windows binaries here.

tangle [switches] inputfile

Input: .node (point set), .poly (planar straight-line graph), or FEMM problem files (.fem, .fee, .feh, .fec).

Output: <base>.1.node, <base>.1.ele, and optionally .edge, .neigh, .poly, .pbc.

# Generate a quality mesh from a PSLG with 30-degree minimum angle
tangle -pq30 myExample

# Same, with area constraint and edge output
tangle -pq30 -a0.1 -e myExample

Algorithms implemented by tangle include the following.

The initial triangulation uses the Bowyer-Watson algorithm (Bowyer 1981, Watson 1981), an incremental insertion method. Points are inserted one at a time into a super-triangle that encloses the entire point set. For each new point, the algorithm finds all existing triangles whose circumcircle contains the point (the “cavity”), removes them, and re-triangulates the resulting polygonal hole with triangles connecting the new point to each edge of the cavity boundary. Points are sorted by x-coordinate before insertion so that a walking point-location strategy (starting from the last inserted triangle) is efficient.

After the initial Delaunay triangulation, PSLG segments that are not already present as edges are enforced to produce a constrained Delaunay triangulation. The algorithm identifies triangles crossed by each missing segment, removes them, and re-triangulates the two polygonal regions on either side using a sweep approach. Segments that cannot be enforced in one pass are split at their midpoint and the process repeats. If a segment still cannot be enforced after repeated splitting, tangle reports an error and writes no output, rather than emit a mesh with a boundary in the wrong place. Constrained edges are then respected as barriers during all subsequent operations.

Holes are removed by Breadth-First Search (BFS) flood-fill from seed points specified in the .poly file, stopping at constrained edges. Region attributes and per-region area constraints are assigned by a similar flood-fill from region seed points.

Bad triangles (minimum angle below the threshold, or area above the constraint) are refined using Ruppert's Delaunay refinement algorithm (Ruppert 1995) with Üngör off-center point placement (Üngör 2004). Instead of inserting at the circumcenter, the off-center method places the new vertex closer to the shortest edge of the bad triangle, producing better-shaped elements and faster convergence.

Bad triangles are processed shortest-edge first. New vertices are inserted into the existing triangulation using Lawson's incremental insertion algorithm (Lawson 1977), which splits the containing triangle and restores the Delaunay property through local edge flips.

Before inserting a new vertex, the algorithm checks whether it would encroach on any constrained segment (i.e., fall inside the segment's diametral circle), or whether the candidate is hidden from its triangle behind a constrained segment — a straight-line visibility walk that also catches circumcenters lying outside the domain entirely. In either case the segment is split at its midpoint instead. This guarantees termination: segment splitting resolves encroachment, and the off-center placement ensures that new vertices do not create progressively worse triangles.

LFS values can be specified per-vertex and per-segment in the .poly file (see polyformat.txt for syntax). During quality refinement, a triangle is considered “bad” not only if its minimum angle or area violates the global thresholds, but also if any of its edges exceeds the LFS value inherited from its endpoints. This provides spatially varying mesh density control: regions near vertices or segments with small LFS values get fine elements, while the rest of the domain stays coarse.

Arc segments define circular arcs between two existing vertices by specifying the subtended angle and a maximum discretization angle. During input processing, each arc is replaced by a chain of straight constrained segments (chord segments) connecting points computed along the arc. Interior vertices created by the discretization are added to the point set and participate normally in the Delaunay triangulation and all subsequent stages.

The arc center is determined from the chord and subtended angle, with the center placed to the left of the chord as traversed from endpoint 1 to endpoint 2. The number of chord segments is ceil(angle / max_seg_angle).

PBC definitions pair two boundary chains identified by their segment boundary markers. During quality refinement, PBC-paired segments are split synchronously: when the refinement algorithm splits a segment on one chain (e.g. because a new Steiner point encroaches on it), it also splits the partner segment on the paired chain at the corresponding parametric position. This maintains a one-to-one node correspondence between the two boundaries throughout refinement.

After meshing, the final node pairs are written to a .pbc output file. This is useful for assembling finite element systems with periodic or anti-periodic coupling between boundaries — for example, exploiting rotational symmetry in electric machines, or using the Kelvin transformation to map unbounded exterior domains onto a finite mesh.

When given a .fem, .fee, .feh, or .fec file, tangle reads the FEMM problem definition directly: node coordinates, line segments, arc segments (with automatic chord discretization), block labels (mapped to regions with area constraints), hole points (<No Mesh> blocks), and boundary properties. SmartMesh and per-segment MaxSideLength constraints are translated to LFS values. PBC and AGE boundary types are detected from BdryFormat codes and handled automatically.

When the -R flag is given (or automatically for FEMM files), tangle reorders the mesh nodes using the Cuthill-McKee algorithm (Cuthill & McKee 1969) to reduce the bandwidth of the resulting finite element stiffness matrix.

The algorithm builds a node adjacency graph from the element topology, then performs a breadth-first traversal starting from the node with minimum connectivity. At each step, unvisited neighbors are appended in order of ascending connectivity (degree), which tends to keep nearby nodes close together in the new numbering. Disconnected components are handled by restarting the BFS from the lowest-connectivity unvisited node.

After renumbering, elements are sorted by the sum of their vertex indices so that elements sharing nodes are stored adjacently. Together, these two reorderings improve cache locality and reduce the bandwidth of sparse matrix storage formats used by direct and iterative FEM solvers.

The geometric predicates orient2d (point-line orientation) and inCircle (point-in-circumcircle) use a floating-point filter with an extended precision arithmetic fallback. Each predicate first evaluates a fast double-precision (~53-bit) computation, then checks the result against a forward error bound derived from the operand magnitudes (Higham 2002). When the result is too close to zero relative to the bound to determine the sign reliably, the computation falls back to extended-precision arithmetic. The extended precision is provided by a custom floating-point type, float256, implemented as one IEEE double (carrying the sign, exponent, and top 52 mantissa bits) plus three unsigned 64-bit integers (uint64_t) that extend the mantissa to ~245 bits.

This two-tier approach (double → float256) avoids the robustness failures that plague double-precision floating-point implementations while keeping the common case fast. Over 99% of predicate calls are resolved at double precision; only the rare near-degenerate cases fall through to the extended precision tier.

Tangle uses the same file formats as Triangle for .node, .poly, .ele, .edge, and .neigh files. See the Triangle documentation for the base format, and polyformat.txt for tangle's extensions (LFS columns, arc segments, PBC definitions, and .pbc output).

Tangle is an independent implementation distributed under the MIT License.

• tangle-0.4.18
• tangle-0.4.17
• tangle-0.4.12
• tangle-0.4.2
• tangle-0.2