Implicit surface mesh creation using Meshes.jl

by

Summary

The core issue is not a lack of capability in Meshes.jl, but a misunderstanding of the mathematical topology of the implicit surface defined by $f(u,v)$. The equation $w = \pm \sqrt{1 – v^2 – (u – v^2)^2}$ describes a surface with a pinch point singularity (specifically, a “pinched torus” or “self-intersecting” topology). It is not a smooth, closed manifold where you can simply “stitch” two halves. Attempting to create a coherent surface mesh directly from this implicit definition without addressing the singularity leads to degenerate triangles or invalid geometry.

Root Cause

The primary cause of failure when trying to stitch the upper and lower halves is the topological singularity at the boundary of the domain. The implicit function is only real-valued where $1 – v^2 – (u – v^2)^2 \ge 0$.

  • The Singularity: At the boundary of this domain, the gradient of the function goes to infinity. The upper and lower halves meet exactly at the boundary curve (where $w=0$).
  • Mathematical Reality: This surface is not a smooth manifold at the “equator.” It behaves like a cone or a pinch.
  • Algorithmic Failure: Standard implicit surface mesher (like Marching Cubes or its variants) expect a smooth gradient. When the gradient blows up or the surface folds back onto itself, the resulting mesh often has zero-area triangles or incorrect normals at the pinch.

Why This Happens in Real Systems

In scientific computing and PDE simulation, geometric validity is non-negotiable. This issue arises frequently when converting implicit mathematical descriptions to explicit computational meshes.

  • Manifold Constraints: PDE solvers (FEM, FVM) require the domain boundary to be a closed manifold (watertight). A mesh containing a pinch point is topologically distinct from a smooth surface; it has singularities where the angle sum around a vertex is not $2\pi$.
  • Mesh Generation Heuristics: Standard mesh generators (Delaunay, Marching Cubes) are designed for “well-behaved” geometry. When the input function has sharp corners or self-intersections (even pinches), these algorithms produce non-conforming or degenerate elements.
  • Discretization Error: Attempting to force a mesh over a singularity introduces massive discretization error. The PDE solver will likely fail to converge or produce garbage results at the pinch point because the basis functions cannot represent the singular flux.

Real-World Impact

  • Solver Instability: The PDE simulation will likely crash due to “negative Jacobian” errors (zero-volume elements) at the pinch point.
  • Physics Artifacts: Even if the solver runs, the singularity acts as a numerical sink or source, creating artificial spikes in the solution variable that do not represent the actual physics.
  • Wasted Development Time: Engineers spend days trying to “fix” the mesh stitching rather than addressing the underlying geometric definition.

Example or Code

The user’s attempt to stitch top and bottom likely fails because the boundaries are topologically curves that share points but do not form a smooth ring. In Meshes.jl, the correct approach to handle implicit surfaces with potential singularities is usually to convert the parametric domain into a parametric mesh and map it, rather than relying on a pure implicit solver that might choke on the pinch.

However, if one must stick to the implicit definition, one usually has to “fatten” the singularity or isolate it. But the cleanest way in Julia is often explicit parametric evaluation.

Here is how you might define the geometry explicitly to avoid the implicit solver’s issues, though you will still face the singularity at the rim:

using Meshes, LinearAlgebra

# Define the implicit function
function f(u, v)
    val = 1 - v^2 - (u - v^2)^2
    if val < 0
        return NaN
    end
    return sqrt(val)
end

# Discretize the parameter space (u, v)
# We must avoid the boundary where val = 0 to prevent NaNs in normals
n = 50
u_range = range(-1.5, 1.5, length=n)
v_range = range(-1.2, 1.2, length=n)

points = Point3{Float64}[]
triangles = []

# Map parameter to 3D space
# Note: We generate a grid and explicitly check connectivity
p_map = Matrix{Union{Point3, Nothing}}(undef, n, n)

for (i, u) in enumerate(u_range)
    for (j, v) in enumerate(v_range)
        w = f(u, v)
        if isnan(w)
            p_map[i, j] = nothing
        else
            # Flip normal if needed for "bottom" half? 
            # The equation sqrt gives positive w. 
            # The implicit is usually S = ±sqrt(...), meaning top and bottom.
            # For a coherent mesh, we need to decide if we want a closed surface 
            # (which this equation doesn't strictly support smoothly) or just the top.
            p_map[i, j] = Point3(u, v, w)
        end
    end
end

# Creating a mesh is tricky here because the domain is simply connected 
# but the surface has a boundary (the pinch). 
# You cannot "stitch" top and bottom because the bottom is just the reflection 
# through the pinch curve.

How Senior Engineers Fix It

Senior engineers address this by recognizing that the mathematical model is the problem, not the meshing code.

  1. Geometry Reformulation: Instead of fighting the pinch, modify the implicit function slightly to create a “smoothed” manifold.
    • Technique: Add a small epsilon: $w = \sqrt{1 – v^2 – (u – v^2)^2 + \epsilon}$. This rounds off the tip, turning the pinch into a smooth tube, making it meshable.
  2. Parametric Mapping: Abandon implicit meshing. Treat the surface as a parametric function $(u, v) \mapsto (x, y, z)$.
    • Action: Create a standard 2D mesh in the $(u, v)$ domain (e.g., a Delaunay triangulation of the valid region). Then map every node to 3D space. This guarantees connectivity without intersection errors.
  3. Domain Truncation: If simulating a PDE, realize that the singularity is a boundary condition. Mesh up to the singularity (shrinking the domain slightly) and apply specific boundary conditions there, rather than trying to mesh through it.

Why Juniors Miss It

  • Over-reliance on Tooling: Juniors often assume that a library like Meshes.jl should “just work” on any mathematical function. They fail to realize that meshing is an approximation of geometry, and if the geometry is topologically broken (singular), the approximation will be too.
  • Visual Confirmation Bias: If they plot the points in 3D, they might look like a surface. They don’t immediately see the zero-volume elements at the pinch that the solver will see.
  • Literal Interpretation: They interpret “stitching top and bottom” as a geometric operation (union) rather than a topological one. They don’t realize that the two halves are conjoined at a curve, not attached along a surface.