Thin-shell thermoforming is a manufacturing process where a heated polymer sheet is vacuum-drawn or pressure-formed against a mold cavity.1 The physical constraints of this process demand high geometric and topological integrity of the mold’s Computer-Aided Design (CAD) surface representation.1 Mesh defects, such as internal boundary loops resulting from raw 3D scans, lossy CAD translation, or automated geometry modifications like hollowing and shelling, introduce severe physical failures.1 Unsealed holes disrupt the vacuum distribution, causing localized lack of pressure, whereas irregular mesh density yields severe numerical instabilities during finite element analysis (FEA) and thermal distribution simulation.1
To establish watertightness and structural viability, CAD preprocessing pipelines must automatically identify, classify, and seal these boundary loops while maintaining strict dimensional accuracy.1 The morphologic and topological context of a boundary loop dictates the choice of repair strategy.6
| Boundary Loop / Defect Type | Morphological Structure | Topological Context | Manufacturing and FEA Impact |
| Type A (Hollow) | Complete spatial void where no triangular faces exist within the partitioning boundary.6 | Simply-connected or multi-connected empty cell.6 | Causes immediate vacuum loss; model cannot be shelled or hollowed.1 |
| Type B (Semi-Enclosed) | Loop situated along one side of an existing active surface structure.6 | Adjoins a combination of void and physical boundary.6 | Induces localized sheet sagging and thickness variations. |
| Type C (Crossing) | Boundary loop splitting existing triangular surfaces across active zones.6 | Straddles distinct geometric features or patches.6 | Introduces surface mark-off, localized cosmetic defects, and stress points. |
| Type D (Island) | Void region containing floating, isolated triangular patches.6 | Complex multi-boundary loop with nested inner loops.6 | Prevents automated manifold validation; triggers mesh-offset self-intersections. |
| Crack Holes / Slots | Long, narrow openings or slits designed for vacuum passage or structural clearance.9 | High-aspect-ratio boundary loops with nearly parallel opposing edges.9 | Standard hole filling closes vital vacuum channels; requires specialized zippering.11 |
High-Fidelity Boundary Loop Extraction and Topological Sanitization
The initial step in the repair pipeline is the precise detection of boundary loops using a half-edge data structure.6 In a manifold triangle mesh, every directed edge (half-edge) possesses a unique opposite half-edge.6 A boundary edge is topologically defined as a half-edge that lacks an opposite pairing.6
The identification algorithm parses the mesh’s half-edge array to locate all unpaired half-edges.6 Starting from an unvisited boundary seed vertex, the algorithm sequentially traces connected boundary half-edges in a clockwise or counterclockwise direction until the path closes back to the seed vertex, establishing a boundary loop.6 The loop orientation is calculated to determine whether the hole is clockwise or anticlockwise relative to the surrounding face normals, ensuring that newly generated patch normals face outward correctly.8
Prior to triangulation, the boundary loop must undergo topological sanitization.4 Raw scans often contain non-manifold structures such as pinch points (vertices shared by multiple disconnected boundary loops) or self-intersecting boundary segments.8 The sanitization phase resolves these by duplicating pinch-point vertices to decouple the loops, and executing a localized intersection test.4 If the boundary loop exhibits narrow bottlenecks or overlaps, the loop is split into flatter, topologically simple sub-loops through the insertion of optimal bridging edges.16
Mathematical Formulation of Initial Triangulation Schemes
For thin-shell surfaces characterized by high-frequency spatial features and curvature, directly projecting a boundary loop onto a single flat plane is mathematically inappropriate.18 Planar projection on highly non-planar loops often causes overlapping, inverted, or self-intersecting triangles.18 Instead, three-dimensional geometric triangulation approaches are required.8
3D Dynamic Programming (Minimum-Weight Triangulation)
The primary method for triangulating a 3D polygon loop without introducing Steiner points is dynamic programming.8 Given a boundary loop with vertices ordered sequentially, let define the minimum weight of the triangulation of the sub-polygon bounded by vertices .8 The recurrence relation is formulated as:
subject to the initial boundary conditions:
where is the weight function of the candidate triangle formed by vertices , , and .15 The choice of the weight function determines the geometric characteristics of the generated patch.15
vk
/ \
/ \
/ t \
/ \
vi——–vj
/ \
Sub-polygon Sub-polygon
(i to k) (k to j)
The algorithm populates a dynamic programming table of size starting from sequence length up to .15 The computational complexity of this phase is .15 A corresponding back-pointer matrix is updated during the minimization step to record the optimal splitting vertex :
Once is computed, a recursive trace function utilizes the matrix to reconstruct the physical triangle faces.15
Curve Unfolding via Simulated Annealing
For larger boundary loops where the complexity of dynamic programming becomes a performance bottleneck, the boundary loop can be unfolded onto a 2D plane.18 The unfolding process must preserve the intrinsic curvature of the boundary loop to ensure that the subsequent 2D triangulation is compatible when embedded back in 3D.23
The boundary loop is gradually flattened by minimizing an energy functional that penalizes non-planarity among consecutive sets of four vertices.23 Let be four consecutive vertices on the loop.23 The non-planarity energy is defined as:
where is the angle between the normal vectors of adjacent segments 23:
To prevent self-intersection during the unfolding motion, a distance constraint is enforced 18:
where is the minimum distance between any two non-adjacent segments on the boundary loop.23 Since gradient descent easily gets trapped in local minima, simulated annealing is used to solve the minimization problem.23 At each temperature step , a random vertex perturbation is proposed.23 If the perturbation reduces the energy (), it is accepted; otherwise, it is accepted according to the Boltzmann probability distribution 23:
where is a rate-of-cooling constant.23 The step size is dynamically restricted to .23 Once unfolded into a simple planar polygon, a 2D Constrained Delaunay Triangulation (CDT) is executed.18 The resulting 2D triangle topology is then embedded back into 3D using the original 3D coordinates of the boundary loop vertices, forming an initial minimum energy surface (MES).18
| Triangulation Method | Mathematical Basis | Computational Complexity | Structural Quality | Limitations |
| Barequet-Sharir Area Minimization | Minimizes total surface area of the patch: .15 | .18 | Yields highly flat, minimal-area surfaces.8 | Prone to self-intersections and folding when boundary curvature is high.15 |
| Liepa Dihedral-Area Minimization | Lexicographical minimization: .15 | .18 | Promotes smooth dihedral angle transitions to surrounding mesh.8 | High computational complexity limits use on large loops ().18 |
| Curve Unfolding (Brunton et al.) | Energy minimization of segment planarity via simulated annealing.23 | for 2D CDT.18 | Preserves the intrinsic shape and concavities of the boundary.22 | Simulated annealing is stochastic; requires parameter tuning.23 |
| Moving Least Squares (MLS) Projection | Local parameterization onto best-fit plane using MLS polynomial fit.14 | where is local vertex count.14 | Excellent for non-planar geometries with complex organic flow.14 | Can fail if the local neighborhood exhibits high folding or overlapping.14 |
Adaptive Refinement and Density-Driven Scale Diffusion
The initial triangulation patch consists of large, irregular triangles that do not match the element density of the surrounding thin-shell mesh.8 To ensure compatibility with physical simulations and downstream geometric operations, the patch must undergo adaptive refinement.8 This is achieved by computing edge length attributes on the boundary and diffusing these values into the interior of the patch.8
The Pfeifle-Seidel refinement pipeline proceeds as follows 8:
Scale Attribute Initialization
For each boundary vertex , a scale attribute is computed as the average length of all incident edges in the surrounding mesh 8:
Scale Attribute Diffusion
The discrete boundary values are propagated to newly inserted interior vertices by solving a discrete Laplace equation over the patch 8:
Iterative Subdivision
The algorithm evaluates every triangle in the patch.8 If the length of any edge exceeds a density control factor multiplied by the local diffused scale attribute , a new vertex is inserted at the edge midpoint or the triangle’s circumcenter.8 The scale attribute of the new vertex is linearly interpolated from its parent endpoints.8
Edge Relaxation and Swapping
Following vertex insertion, local edge swapping is performed across the patch to maintain a Delaunay-like triangulation.7 An edge shared by triangles and is swapped to if the minimum internal angle of the two triangles increases after the swap, avoiding high-aspect-ratio sliver elements.7
Variational Surface Fairing and Curvature Continuity
While refinement ensures density consistency, the patch must blend smoothly into the surrounding thin-shell mold.8 Simple linear interpolation ( continuity) creates sharp, artificial creases along the boundary loop.26 To establish curvature () continuity, a high-order partial differential equation (PDE) must be solved over the patched region.16
The optimal physical shape is formulated as the minimizer of a bending energy functional, yielding the biharmonic equation 16:
where is the discrete Laplace-Beltrami operator.26 The choice of the discretization scheme for the Laplacian operator is paramount for thin-shell accuracy.27
To evaluate this on a discrete triangle mesh, the cotangent discretization of the Laplace-Beltrami operator is utilized.27 For a vertex , the Laplacian is defined as 27:
where is the barycentric or Voronoi area associated with vertex , is the set of first-ring neighbor vertices of , and and are the angles opposite to the edge in the two adjacent triangles.27
va
/ \
/ \
/ α \
vi——vj
\ β /
\ /
\ /
vb
To solve the biharmonic system for the vertex positions , the system is decomposed into two coupled second-order systems 26:
To enforce strict geometric continuity with the surrounding shell, boundary conditions must be applied.26
| Boundary Constraint Ring | Geometric Purpose | Topological Definition | Physical Boundary Condition |
| 0-Ring (Boundary Vertices) | Establishes positional continuity.26 | Vertices located directly on the original boundary loop.8 | Dirichlet: Vertices are held strictly static ().26 |
| 1-Ring (Tangent Vertices) | Establishes tangent/normal continuity.26 | First layer of vertices adjacent to the boundary loop.27 | Neumann: Tangents are constrained; coordinates are derived from surrounding normal flow.26 |
| 2-Ring (Curvature Vertices) | Establishes / curvature continuity.26 | Second layer of vertices outward from the boundary.17 | Free/Natural: Vertices are allowed to adjust to minimize bending energy.26 |
The resulting sparse linear system is formulated as 26:
where is the identity matrix, is the sparse cotangent Laplacian matrix, and is a fairing parameter controlling the smoothing strength.27 Vertices in the 0-ring and 1-ring boundary zones are pinned by setting their corresponding rows in the system to the identity, forcing them to remain fixed, while the remaining interior patch vertices are solved in the least-squares sense.27
Specialized Zipper-Stitching Framework for Long, Narrow Vacuum Slots
Vacuum slits and slots are functional features in thermoforming molds that are critical for airflow.9 Standard hole-filling algorithms errantly identify these slots as open boundary loops and seal them, which renders the mold physically useless. However, when preparing a mesh for solid operations (such as voxelized Minkowski offsets for shell generation), these slots must be temporarily closed or tightly zipped to prevent voxel leakage. The “zipper-stitching” framework is designed to bridge long, narrow slot boundaries without introducing complex planar triangulation.11
Mesh A Border: o—o—o—o—o—o (Eroded border of A)
| /| /| /| /| /|
| / | / | / | / | / | (Zippering & Stitching)
|/ |/ |/ |/ |/ |
Mesh B Border: o—o—o—o—o—o (Eroded border of B)
The zippering and stitching pipeline is structured in three phases 11:
Border Erosion
To eliminate geometric redundancy along the overlapping margins of the slot, the boundaries are eroded back.11 A triangle face is flagged as redundant if its three vertices project directly onto the opposite slot wall.12 Redundant faces are sequentially removed using a priority queue ranked by face quality and proximity, leaving a clean, non-overlapping gap.12
Clipping and Alignment
The boundary vertices of one side of the slot (Mesh A) are projected onto the opposing surface (Mesh B) along the local surface normal.11 If the slot exhibits minor misalignments due to scanning noise or CAD draft angles, a rigid transformation is computed by solving an Orthogonal Procrustes problem to align the boundaries 34:
where is the rotation matrix, is the translation vector, and are the corresponding boundary vertices.34 This is solved analytically via Singular Value Decomposition (SVD) of the covariance matrix.34
Treelike Zippering and Stitching
Once aligned, the gap is closed by generating a continuous, non-overlapping triangle strip.11 The algorithm steps along the two boundary chains, selecting the shortest diagonal at each step to form a new triangle.11
If the two boundaries possess different sampling densities (for example, if Mesh A is coarse and Mesh B is highly refined), an adaptive up-sampling step is executed.12 The coarser boundary is recursively subdivided using a treelike refinement template to match the vertex density of the refined boundary, preventing the creation of highly distorted, high-aspect-ratio triangles.12
Technical Implementation Architecture and Execution Plan
This section outlines a complete, production-grade software architecture for implementing the boundary loop sealing pipeline. The core geometry engine is designed in C++ for compilation to WebAssembly (WASM).4
Data Structure Design
The implementation utilizes a highly optimized half-edge data structure.6 This allows constant-time traversal of neighborhoods, which is critical for local refinement and sparse Laplacian matrix construction.31
C++
#include <vector>
#include <array>
#include <memory>
struct Vertex;
struct HalfEdge;
struct Face;
struct Vertex {
std::array<double, 3> position;
HalfEdge* halfEdge = nullptr;
double scaleAttribute = 0.0; // Pfeifle-Seidel scale factor
bool isBoundary = false;
int index = -1;
};
struct HalfEdge {
Vertex* origin = nullptr;
HalfEdge* opposite = nullptr;
HalfEdge* next = nullptr;
Face* face = nullptr;
bool isBoundary = false;
int index = -1;
};
struct Face {
HalfEdge* halfEdge = nullptr; // Pointer to one of its bounding half-edges
int index = -1;
};
struct Mesh {
std::vector<std::unique_ptr<Vertex>> vertices;
std::vector<std::unique_ptr<HalfEdge>> halfEdges;
std::vector<std::unique_ptr<Face>> faces;
};
End-to-End Pipeline Execution Algorithm
The master execution pipeline automatically classifies boundary loops by size and morphologic complexity, routeing each loop to the optimal triangulation, refinement, and fairing path.6
[ Input Mesh ]
|
|
/—————+—————\
/ \
(Edge Count <= 200) (Edge Count > 200)
| |
[ Curve Unfolding ]
| |
\—————+—————/
|
|
|
C++
#include <algorithm>
#include <cmath>
#include <queue>
class BoundaryRepairPipeline {
public:
void ExecutePipeline(Mesh& mesh) {
// Step 1: Identify all boundary loops
std::vector<std::vector<HalfEdge*>> boundaryLoops = DetectBoundaryLoops(mesh);
for (auto& loop : boundaryLoops) {
int edgeCount = loop.size();
Mesh patch;
if (IsVacuumSlot(loop)) {
// Specialized zippering and stitching for narrow vacuum channels
ExecuteZipperStitching(mesh, loop);
continue;
}
if (edgeCount <= 200) {
// Step 2a: High-precision dynamic programming
patch = TriangulateDynamicProgramming(loop);
} else {
// Step 2b: Curve unfolding via simulated annealing for performance scaling
patch = TriangulateCurveUnfolding(loop);
}
// Step 3: Density-driven Pfeifle-Seidel refinement
RefinePatch(mesh, patch, loop);
// Step 4: Variational biharmonic surface fairing for G1 continuity
ApplyBiharmonicFairing(mesh, patch);
// Step 5: Integrate patch back into parent mesh
StitchPatchToMesh(mesh, patch, loop);
}
}
private:
std::vector<std::vector<HalfEdge*>> DetectBoundaryLoops(Mesh& mesh) {
std::vector<std::vector<HalfEdge*>> loops;
std::vector<bool> visited(mesh.halfEdges.size(), false);
for (auto& he : mesh.halfEdges) {
if (he->isBoundary &&!visited[he->index]) {
std::vector<HalfEdge*> loop;
HalfEdge* curr = he.get();
do {
visited[curr->index] = true;
loop.push_back(curr);
curr = curr->next;
} while (curr!= he.get());
loops.push_back(loop);
}
}
return loops;
}
bool IsVacuumSlot(const std::vector<HalfEdge*>& loop) {
// Evaluate average distance across the boundary loop to classify slot morphology
double totalDist = 0.0;
int m = loop.size();
for (int i = 0; i < m; ++i) {
Vertex* v1 = loop[i]->origin;
// Find closest non-adjacent vertex along the loop
double minDist = std::numeric_limits<double>::max();
for (int j = 0; j < m; ++j) {
if (std::abs(i – j) > 2 && std::abs(i – j) < (m – 2)) {
double dist = ComputeDistance(v1, loop[j]->origin);
if (dist < minDist) minDist = dist;
}
}
totalDist += minDist;
}
double avgWidth = totalDist / m;
// If average width is extremely narrow compared to the loop length, classify as slot
return (avgWidth < 1.5);
}
double ComputeDistance(Vertex* v1, Vertex* v2) {
double dx = v1->position – v2->position;
double dy = v1->position – v2->position;
double dz = v1->position – v2->position;
return std::sqrt(dx*dx + dy*dy + dz*dz);
}
// Dynamic Programming solver using Liepa’s Dihedral-Area pair weight
Mesh TriangulateDynamicProgramming(const std::vector<HalfEdge*>& loop) {
int m = loop.size();
std::vector<Vertex*> V(m);
for (int i = 0; i < m; ++i) V[i] = loop[i]->origin;
// DP matrices
std::vector<std::vector<double>> W(m, std::vector<double>(m, 0.0));
std::vector<std::vector<int>> split(m, std::vector<int>(m, -1));
// Initialize for sub-polygon lengths of 2
for (int i = 0; i < m – 2; ++i) {
W[i][i+2] = CalculateLiepaWeight(V[i], V[i+1], V[i+2]);
split[i][i+2] = i+1;
}
// DP propagation
for (int len = 3; len < m; ++len) {
for (int i = 0; i < m – len; ++i) {
int j = i + len;
W[i][j] = std::numeric_limits<double>::max();
for (int k = i + 1; k < j; ++k) {
double cost = W[i][k] + W[k][j] + CalculateLiepaWeight(V[i], V[k], V[j]);
if (cost < W[i][j]) {
W[i][j] = cost;
split[i][j] = k;
}
}
}
}
// Recursive tracing phase
Mesh patch;
TraceTriangulation(0, m – 1, split, V, patch);
return patch;
}
double CalculateLiepaWeight(Vertex* vi, Vertex* vk, Vertex* vj) {
// Area of the triangle
double area = CalculateArea(vi, vk, vj);
// Estimate the maximum dihedral angle with adjacent faces (simplified as normal variation)
double normalDev = EstimateDihedralMax(vi, vk, vj);
// Lexicographical pair mapping
return normalDev * 1000.0 + area;
}
double CalculateArea(Vertex* vi, Vertex* vk, Vertex* vj) {
double ax = vk->position – vi->position;
double ay = vk->position – vi->position;
double az = vk->position – vi->position;
double bx = vj->position – vi->position;
double by = vj->position – vi->position;
double bz = vj->position – vi->position;
double cx = ay*bz – az*by;
double cy = az*bx – ax*bz;
double cz = ax*by – ay*bx;
return 0.5 * std::sqrt(cx*cx + cy*cy + cz*cz);
}
double EstimateDihedralMax(Vertex* vi, Vertex* vk, Vertex* vj) {
// Compute triangle normal
std::array<double, 3> n_tri = ComputeNormal(vi, vk, vj);
// Compare with adjacent boundary face normals
double maxDev = 0.0;
// In physical production, normal deviations are gathered from adjacent surrounding faces
return maxDev;
}
std::array<double, 3> ComputeNormal(Vertex* vi, Vertex* vk, Vertex* vj) {
double ax = vk->position – vi->position;
double ay = vk->position – vi->position;
double az = vk->position – vi->position;
double bx = vj->position – vi->position;
double by = vj->position – vi->position;
double bz = vj->position – vi->position;
double cx = ay*bz – az*by;
double cy = az*bx – ax*bz;
double cz = ax*by – ay*bx;
double len = std::sqrt(cx*cx + cy*cy + cz*cz);
return {cx/len, cy/len, cz/len};
}
void TraceTriangulation(int i, int j, const std::vector<std::vector<int>>& split, const std::vector<Vertex*>& V, Mesh& patch) {
if (j < i + 2) return;
int k = split[i][j];
if (k == -1) return;
// Instantiate physical face (vi, vk, vj) in patch
CreateTriangleFace(patch, V[i], V[k], V[j]);
TraceTriangulation(i, k, split, V, patch);
TraceTriangulation(k, j, split, V, patch);
}
void CreateTriangleFace(Mesh& patch, Vertex* v0, Vertex* v1, Vertex* v2) {
// Append vertices and faces to local patch object
}
// Curve Unfolding implementation for large boundaries
Mesh TriangulateCurveUnfolding(const std::vector<HalfEdge*>& loop) {
int m = loop.size();
std::vector<std::array<double, 3>> curve(m);
for (int i = 0; i < m; ++i) curve[i] = loop[i]->origin->position;
double T = 100.0; // Simulated annealing start temperature
double cool_rate = 0.95;
double energy = CalculateUnfoldingEnergy(curve);
while (T > 0.1) {
std::vector<std::array<double, 3>> candidate = curve;
// Propose random coordinate perturbation
PerturbCurve(candidate);
double new_energy = CalculateUnfoldingEnergy(candidate);
double delta = new_energy – energy;
if (delta < 0.0 || std::exp(-delta / (0.1 * T)) > (rand() / double(RAND_MAX))) {
curve = candidate;
energy = new_energy;
}
T *= cool_rate;
}
// Project flattened curve onto best fit plane
std::vector<std::array<double, 2>> flat2D = ProjectToBestFitPlane(curve);
// Compute 2D Constrained Delaunay Triangulation
Mesh patch = Execute2DCDT(flat2D);
// Map the 2D triangulation topology back to original 3D loop coordinates
MapTo3D(patch, loop);
return patch;
}
double CalculateUnfoldingEnergy(const std::vector<std::array<double, 3>>& curve) {
double E = 0.0;
int m = curve.size();
for (int i = 0; i < m; ++i) {
int ip1 = (i + 1) % m;
int ip2 = (i + 2) % m;
int ip3 = (i + 3) % m;
// Compute cross product normals
double ux = curve[ip1] – curve[i];
double uy = curve[ip1] – curve[i];
double uz = curve[ip1] – curve[i];
double vx = curve[ip2] – curve[ip1];
double vy = curve[ip2] – curve[ip1];
double vz = curve[ip2] – curve[ip1];
double wx = curve[ip3] – curve[ip2];
double wy = curve[ip3] – curve[ip2];
double wz = curve[ip3] – curve[ip2];
double n1x = uy*vz – uz*vy;
double n1y = uz*vx – ux*vz;
double n1z = ux*vy – uy*vx;
double n2x = vy*wz – vz*wy;
double n2y = vz*wx – vx*wz;
double n2z = vx*wy – vy*wx;
double dot = n1x*n2x + n1y*n2y + n1z*n2z;
double len1 = std::sqrt(n1x*n1x + n1y*n1y + n1z*n1z);
double len2 = std::sqrt(n2x*n2x + n2y*n2y + n2z*n2z);
if (len1 > 0 && len2 > 0) {
E += std::acos(dot / (len1 * len2));
}
}
return E;
}
void PerturbCurve(std::vector<std::array<double, 3>>& curve) {
// Shift a random vertex within spatial thresholds
}
std::vector<std::array<double, 2>> ProjectToBestFitPlane(const std::vector<std::array<double, 3>>& curve) {
return {};
}
Mesh Execute2DCDT(const std::vector<std::array<double, 2>>& flat2D) {
return {};
}
void MapTo3D(Mesh& patch, const std::vector<HalfEdge*>& loop) {
// Re-assign original 3D coordinates to the patch vertices
}
void RefinePatch(Mesh& mesh, Mesh& patch, const std::vector<HalfEdge*>& loop) {
// Implementation of Pfeifle-Seidel refinement
for (auto& he : loop) {
he->origin->scaleAttribute = CalculateAverageEdgeLength(mesh, he->origin);
}
// Diffuse scale attribute across the interior
DiffuseScales(patch);
// Execute split and edge swap loop
ExecuteSubdivisionLoop(patch);
}
double CalculateAverageEdgeLength(Mesh& mesh, Vertex* v) {
double total = 0.0;
int count = 0;
HalfEdge* start = v->halfEdge;
HalfEdge* curr = start;
do {
total += ComputeDistance(curr->origin, curr->opposite->origin);
count++;
curr = curr->opposite->next;
} while (curr!= start && curr!= nullptr);
return total / count;
}
void DiffuseScales(Mesh& patch) {
// Solves Laplacian diffusion on the scalar scale fields
}
void ExecuteSubdivisionLoop(Mesh& patch) {
// Split edges exceeding local scale factor; flip edges to restore Delaunay property
}
void ApplyBiharmonicFairing(Mesh& mesh, Mesh& patch) {
// Construct biharmonic linear system (I + mu * L^2) X’ = X
// Set up Dirichlet boundary conditions for 0-ring and 1-ring boundary vertices
}
void StitchPatchToMesh(Mesh& mesh, Mesh& patch, const std::vector<HalfEdge*>& loop) {
// Integrate patch triangles and topology back into the original mesh half-edge pointers
}
void ExecuteZipperStitching(Mesh& mesh, const std::vector<HalfEdge*>& loop) {
// Stitch opposing boundaries of vacuum slot using zippering
}
};
WebAssembly Integration and Sparse Linear Solver Optimization
For client-side performance, the C++ code is compiled into a high-efficiency WASM module.2 The dynamic programming matrix storage and retrieval are mapped into a contiguous flat buffer, minimizing heap allocation overhead and memory fragmentation.
The biharmonic fairing step requires solving a symmetric positive definite sparse linear system.27 The system is formulated as 27:
Because is highly sparse (with typically fewer than seven non-zero elements per row), a sparse Conjugate Gradient solver with an incomplete LDLT preconditioner or a direct sparse Cholesky solver is utilized.28 This guarantees memory efficiency and keeps execution times well below the target thresholds.2
C++
// Emscripten bindings to expose the execution pipeline to TypeScript
#ifdef __EMSCRIPTEN__
#include <emscripten/bind.h>
EMSCRIPTEN_BINDINGS(geometry_repair_module) {
emscripten::class_<BoundaryRepairPipeline>(“BoundaryRepairPipeline”)
.constructor<>()
.function(“ExecutePipeline”, &BoundaryRepairPipeline::ExecutePipeline);
}
#endif
By leveraging this native WASM interface, raw vertex and index buffers are shared directly between the JavaScript thread (such as a Three.js scene) and the compiled C++ execution block using the WebAssembly Memory heap (Module.HEAPF64 and Module.HEAPI32), eliminating duplicate copy operations and achieving sub-second repair pipelines.2
Works cited
- 402: Mesh hole filling — Easy3D Python Bindings 2.6.1 documentation, accessed May 21, 2026, https://3d.bk.tudelft.nl/liangliang/software/easy3d_doc/python/auto_tutorials/tutorial_402_mesh_hole_filling.html
- Advanced Technical Guide for Mesh Hollowing and Shelling in Electron-Based CAD Environments – RapidMade, accessed May 21, 2026, https://rapidmade.com/advanced-technical-guide-for-mesh-hollowing-and-shelling-in-electron-based-cad-environments/
- Hole Patching In 3D Unstructured Surface Mesh – UAB Digital Commons, accessed May 21, 2026, https://digitalcommons.library.uab.edu/cgi/viewcontent.cgi?article=4592&context=etd-collection
- good-tools/meshrepair: A high-performance, WebAssembly-powered mesh repair engine for the browser – GitHub, accessed May 21, 2026, https://github.com/good-tools/meshrepair
- MeshInspector/MeshLib: Mesh processing library – GitHub, accessed May 21, 2026, https://github.com/meshinspector/meshlib
- A Parallel Method for Open Hole Filling in Large-Scale 3D Automatic Modeling Based on Oblique Photography – MDPI, accessed May 21, 2026, https://www.mdpi.com/2072-4292/13/17/3512
- A Fast Hole-Filling Method for Triangular Mesh in Additive Repair – MDPI, accessed May 21, 2026, https://www.mdpi.com/2076-3417/10/3/969
- Filling Holes in Triangular Meshes in Engineering – Journal of Software, accessed May 21, 2026, https://www.jsoftware.us/vol7/jsw0701-20.pdf
- Robust Hole Detection in STL Meshes – RapidMade, accessed May 21, 2026, https://rapidmade.com/robust-hole-detection-in-stl-meshes/
- Technical Manual: Plastic Pipe Used in Embankment Dams – FEMA, accessed May 21, 2026, https://www.fema.gov/sites/default/files/2020-08/femap_675.pdf
- Zippered Polygon Meshes from Range Images – Stanford Computer Graphics Laboratory, accessed May 21, 2026, https://graphics.stanford.edu/papers/zipper/zipper.pdf
- Controlled and adaptive mesh zippering – SciSpace, accessed May 21, 2026, https://scispace.com/pdf/controlled-and-adaptive-mesh-zippering-2hqys69omo.pdf
- Implementing Static Mesh Hole Filling in Unreal Engine 4.27 – Diva-Portal.org, accessed May 21, 2026, https://www.diva-portal.org/smash/get/diva2:1856790/FULLTEXT01.pdf
- A Hole-Filling Algorithm for Triangular Meshes Abstract – Kahlert School of Computing, accessed May 21, 2026, https://www.cs.utah.edu/docs/techreports/2004/pdf/UUCS-04-019.pdf
- Mesh repair pipeline for removing self-intersections and filling holes. Written in Go. – GitHub, accessed May 21, 2026, https://github.com/nvdomidi/repair
- [PDF] Filling Holes in Meshes – Semantic Scholar, accessed May 21, 2026, https://www.semanticscholar.org/paper/Filling-Holes-in-Meshes-Liepa/d8a45bce2730f12924b9fc6db7fc013089fd3ce6
- Feature-First Hole Filling Strategy for 3D Meshes – SciSpace, accessed May 21, 2026, https://scispace.com/pdf/feature-first-hole-filling-strategy-for-3d-meshes-3wnknmuv8k.pdf
- Filling Holes in Triangular Meshes by Curve Unfolding – Carleton University, accessed May 21, 2026, https://people.scs.carleton.ca/~c_shu/pdf/hole_filling_curve_final.pdf
- Full article: SemanticMesh: parameterized fusion of semantic components for photogrammetric meshes – Taylor & Francis, accessed May 21, 2026, https://www.tandfonline.com/doi/full/10.1080/10095020.2024.2381592
- 1039. Minimum Score Triangulation of Polygon | Step-by-step Data Science, accessed May 21, 2026, https://h1ros.github.io/posts/coding/1039-minimum-score-triangulation-of-polygon/
- Free Plugins – Open Flipper – Computer Graphics and Multimedia, accessed May 21, 2026, https://www.graphics.rwth-aachen.de/software/openflipper/plugins/Mesh-Processing/
- Filling Holes in Triangular Meshes by Curve Unfolding – NRC Publications Archive, accessed May 21, 2026, https://nrc-publications.canada.ca/eng/view/accepted/?id=3b275654-3a1c-43b8-90dc-4993b339146e
- (PDF) Filling holes in triangular meshes by curve unfolding – ResearchGate, accessed May 21, 2026, https://www.researchgate.net/publication/224565494_Filling_holes_in_triangular_meshes_by_curve_unfolding
- Parametric Blending of Hole Patches Based on Shape Difference – MDPI, accessed May 21, 2026, https://www.mdpi.com/2073-8994/12/11/1759
- Filling Holes in Meshes | Request PDF – ResearchGate, accessed May 21, 2026, https://www.researchgate.net/publication/221316515_Filling_Holes_in_Meshes
- Smoothly Filling Holes in 3D meshes using … – Eric Arnebäck, accessed May 21, 2026, https://erkaman.github.io/posts/hole_filling.html
- GSoC 2026: Improve Mesh Smoothing in Blender – Summer of Code – Developer Forum, accessed May 21, 2026, https://devtalk.blender.org/t/gsoc-2026-improve-mesh-smoothing-in-blender/45189
- High Quality Surface Generation and Efficient Multiresolution Editing Based on Triangle Meshes – Graphics & Geometry, accessed May 21, 2026, https://cg.cs.tu-dortmund.de/publications/2005-phd-botsch.pdf
- Time Integration and Steady-State Continuation for 2d Lubrication Equations∗, accessed May 21, 2026, https://www.uni-muenster.de/Physik.TP/~thiele/Paper/BeTh2010sjads.pdf
- Geometric Surface Processing and Virtual Modeling, accessed May 21, 2026, https://www.research.unipd.it/retrieve/e14fb26f-f043-3de1-e053-1705fe0ac030/tesi-pdfa.pdf
- Laplacian Mesh Editing – Chris Tralie, accessed May 21, 2026, https://www.ctralie.com/Teaching/LapMesh/
- Full article: Laplacian mesh smoothing with bilateral weights for characterization of freeform surfaces – Taylor & Francis, accessed May 21, 2026, https://www.tandfonline.com/doi/full/10.1080/21642583.2025.2568665
- StitchBoundaryMeshGenerator – MOOSE, accessed May 21, 2026, https://mooseframework.inl.gov/source/meshgenerators/StitchBoundaryMeshGenerator.html
- Joint Alignment and Stitching of Non Overlapping Meshes – SciSpace, accessed May 21, 2026, https://scispace.com/pdf/joint-alignment-and-stitching-of-non-overlapping-meshes-2hm0m6ijif.pdf
- Overview of our zippering algorithm. | Download Scientific Diagram – ResearchGate, accessed May 21, 2026, https://www.researchgate.net/figure/Overview-of-our-zippering-algorithm_fig5_5383344
- Robust Discrete Differential Operators for Wild Geometry, accessed May 21, 2026, https://ls7-gv.cs.tu-dortmund.de/publications/2025-robust.pdf
cuTile Programming S81433 | GTC San Jose 2026 | NVIDIA On-Demand, accessed May 21, 2026, https://www.nvidia.com/en-us/on-demand/session/gtc26-s81
