General concepts » Fitting a mesh interface with SWIFT

Setup and ownership

Problem fits marked mesh facets to a supplied zero level set. It writes accumulated displacement into the constructor's trial-function solution. It does not move vertices, replace mesh transformations, or remesh. The mesh is obtained from the displacement finite element space.

The mesh, finite element space, trial and test functions must outlive the solver. Use matching vector-valued trial and test spaces on the local Eigen backend. Mark the interface facets and prepare the mesh connectivity before solving. The target gradient must be the derivative of the target level set at the displaced points; no target Hessian is required.

MinSTCut is the canonical classifier used by the reconstruction examples. It balances cell-average phase fidelity against the measure of separating facets. Apply its inside and outside sets to cell attributes and its cut set to interface facets before fitting. Classification remains an explicit preparation step, not part of the SWIFT optimization or a boundary condition.

// mesh already has its interface facets marked with Interface.
P1<Math::SpatialVector<Real>, Geometry::LocalMesh> space(mesh, 2);
// u owns the accumulated displacement; v tests increments on the same space.
TrialFunction u(space);
TestFunction v(space);
Adaptation::SWIFT::Problem fitting(u, v);
Adaptation::SWIFT::Parameters parameters;
// The reference spacing stays fixed as the displacement changes.
parameters.model.h = referenceSpacing;
fitting.setParameters(parameters).setInterfaceAttribute(Interface);
// Zero increments keep this boundary at its current position.
fitting += DirichletBC(u, VectorFunction(Real(0), Real(0))).on(FixedBoundary);

// Fit to a vertical line using its exact gradient; no Hessian is needed.
RealFunction phi([](const Geometry::Point& point) {
  return point.x() - Real(0.55);
});
VectorFunction gradient(Real(1), Real(0));
const auto report = fitting.solve(phi, gradient);
// Check report.reason, geometricSup, minJ and maxQRel before using the result.
// Applying displacement to mesh geometry is a separate application step.
const auto& displacement = u.getSolution();

Homogeneous boundary conditions hold the current boundary position fixed: they constrain each increment, not the accumulated displacement to zero. Nonzero increment values and identification conditions are rejected. For another solve on the same space, retain the solver and update its parameters or target. The existing displacement is the next initial guess. Create a new solver when replacing the mesh or finite element space.

Applying the result is a separate application decision. Vertex displacement is appropriate for a linear geometric mesh. Higher-order geometry requires updating its transformations, not merely moving its corner vertices.

For an affine simplicial mesh, Adapt owns the vector P1 space and fitting problem and applies valid displacements in place:

// Prepare connectivity and mark interface facets before constructing adapt.
Adaptation::SWIFT::Adapt adapt(mesh);
Adaptation::SWIFT::Parameters controls;
controls.model.h = referenceSpacing; // Keep the background scale fixed.
adapt.setParameters(controls).setInterfaceAttribute(Interface);
// Use the owned trial function for additional boundary conditions.
auto& trial = adapt.getTrialFunction();
adapt.getProblem() +=
  DirichletBC(trial, VectorFunction(Real(0), Real(0))).on(FixedBoundary);
const auto adapted = adapt.execute(phi, gradient); // Mesh is updated here.

Each execute call resets displacement to zero on the current mesh. It applies a result when qualityBudgetSatisfied is true, including valid best-effort exits, and leaves the mesh untouched otherwise. The report describes the fit relative to that call's starting mesh. The mesh must outlive the workflow; replacing connectivity requires a new workflow. Classification and remeshing remain application responsibilities. getProblem() provides the same metric, boundary-condition and monitor interfaces as an explicitly constructed problem.

Backend selection follows storage and mesh context, not a runtime conversion. An MPI mesh will use PETSc when that specialization is implemented; it never falls back to Eigen. Renaming the API does not provide distributed assembly.

Mathematical model

The displacement metric is : fitting observes normal motion on the marked interface, while distribution couples that motion to the volume. The metric defines the quadratic direction problem; it is not an additional energy minimized over the accumulated deformation. Quality is recovered through affine squared hinges on the increment and nonlinear checks on the accepted mesh, not through a separate shape-curvature term in the metric.

For the current map , set , and . The two bilinear forms are

Here and are the reference interface and volume; is a fixed sampled target-gradient normalization. Fitting uses squared-residual curvature without the level-set Hessian or robust weights. The force and outer Armijo merit instead come from the fixed-scale Welsch objective

The reference measure and robust scale are held fixed during a solve. Thus the quadratic metric preconditions the robust descent; it does not replace the robust objective with a least-squares objective.

Directional Newton scales the predictor before the affine hinge solve: and . The inner model is

where is the squared-hinge penalty on frozen Jacobian and distortion slacks. The penalty and actual-quality checks use the same cached closed-form reference covering. The generating lattice is shifted and contracted without changing its point count; vertices are not retained as supplemental points. A known reference covering radius does not certify continuous mesh quality. Separate positive Jacobian and distortion weights mix equal mass with normalized nonlinear guard penetration at the current and full-predictor geometries. Each family preserves mapped cell mass, and its weights stay fixed during the inner solve. The equal fraction is one half and dimensionless risks are capped at 100. Inverted predictor points receive maximal distortion risk without evaluating distortion. Validation checks every witness without weights. Equal reference weights, normalized to reference volume and mapped to each cell, determine cell mass, not individual adaptive penalty weights. Coverings are not degree-exact quadrature rules and do not certify nonlinear quality between witnesses. Its stationarity residual controls inner convergence. The outer Armijo search then decreases while enforcing the actual sampled quality budget. See the SWIFT class reference for the slack definitions, penalty scaling, Newton system, line search and geometric stopping equations.

For current-coordinate fields, ; the overbars are averages over the entire current domain . They are not per-element means. Both distribution coefficients are independent.

What the solver does

Each outer iteration has four stages:

1. Assemble the robust fitting force and the metric . Fitting makes motion responsive to the target normal. Distribution couples neighboring elements through variation of current strain about its global mean. Coherent affine motions are free, including shear and anisotropic stretching. Spatial strain variation is resisted. This does not permit independent incompatible element motions; shared DOFs transmit motion.
2. Solve for an unconstrained predictor. Directional Newton rescales the physical predictor and its quadratic metric before applying the hinges.
3. Use native NewtonSolver to solve the frozen quadratic model with affine quadratic quality hinges. Inactive hinges need no Newton correction. Inner convergence uses the stationarity residual.
4. Backtrack with Armijo on the robust fitting energy. Every trial must also satisfy the actual sampled Jacobian and distortion limits. Accept the displacement only after both checks pass.

The retained linear adapter supports CG, SparseLU and optional MUMPS. MUMPS retains symbolic analysis when the sparsity pattern is unchanged. Numerical factors are reused only when the assembled operator is unchanged. Global centering is represented by a low-rank subtraction from the sparse core. Direct solvers factor the equivalent augmented system; CG applies the centered operator without a dense matrix. No gauge, mass penalty or stabilization is added. Singular systems remain subject to the selected backend. With both centered weights positive, distribution is coercive in H1 modulo global affine motions on a fixed connected bounded Lipschitz domain. The mesh-scale factor and changing domain preclude an unconditional uniform bound.

Main controls

Model coefficients are grouped under parameters.model. Accuracy settings and work budgets are independent groups under parameters.convergence: tolerance and iterations. Globalization, linear-solver selection and quadrature have their own groups. All values in the following table are defaults except the required reference size.

MUMPS is the default when compiled in; otherwise SparseLU is selected. Set linear.solver and linear.threads explicitly for reproducible runs. The automatic sampling policies are:

The affine policy is calibrated on the tested cases. Higher-order policies are provisional and require refinement checks, particularly when hinges are active. Orders grow for degrees beyond this table; those floors are heuristic policies, not exactness guarantees for nonlinear integrands. Selection is per entity. Nonsimplicial volume elements use the higher-order volume and quality policies even at degree one; surface integration uses the facet's geometry and trace degree. Quadrature orders are selected automatically unless overridden by quadrature.order (example flag --quadrature-order). A known fitting tensor on an affine entity uses the complete integrand's native order propagation. A globally constant target gradient is unchanged by point location and yields order for degree- displacement functions. Elementwise target orders alone cannot be propagated through deformation: a moved facet may cross target cells, and curved surface measures may be nonpolynomial. Such tensors retain the nonpolynomial interface policy, as do the Welsch force and energy; these two use the same integration rule. quadrature.validation (--quadrature-validation) independently controls geometric-response sampling and is not reduced by --quadrature-order. The surface and volume settings use --quadrature-surface and --quadrature-volume. Independent quality checks use --sampling-subdivision, interpreted as generating lattice subdivisions, not polynomial exactness. Tensor cells use Cartesian grids, wedges use triangular grids times segments, and pyramids use shrinking square layers. Default triangle/tetrahedron counts are for affine and for . Cell vertices are always checked, and common integration overrides do not change independent quality sampling. These checks remain sampling, not an everywhere certificate on curved cells. Example flags follow the C++ parameter hierarchy: dots become hyphens and camel-case words are separated by hyphens. Thus sampling.subdivision maps to --sampling-subdivision. Model flags include --model-fit, --model-distribution-deviatoric, --model-distribution-divergence, --model-hinge, --model-jacobian and --model-distortion. Iteration flags are --convergence-iterations-outer and --convergence-iterations-inner. Linear settings use --linear-solver, --linear-threads and --convergence-tolerance-linear-relative. Legacy spellings are rejected.

Reading results

Use report.reason rather than parsing a message. A geometric target hit is SWIFT::Report::Reason::GeometricTarget. Budget, small-step and small-energy exits are best-effort results, not geometric convergence. The stable text token is available through report.getReasonString().

The main geometric response is geometricSup, the maximum sampled normalized level-set residual on the complete interface, including its vertices. Normalization divides the absolute level-set value by the local target-gradient magnitude, giving a local distance estimate rather than a raw level-set value. The automatic target is , using the reference spacing and the interface finite element order. This is not a certified Hausdorff error; a single target hit does not establish a refinement order.

The report's geometricConstant is the dimensionless sampled score

The background reference size remains fixed throughout fitting. This normalization is independent of a user-specified stopping tolerance. Geometry traces and canonical example responses expose the score as geom_c.

Monitor energy, geometricSup, geometricConstant, geometricSupTarget, minJ, maxQRel, iterations, innerIterations, lastInnerIterations, maxInnerIterations and innerRelativeResidual. The report's innerIterations is the total correction count, whereas parameters.convergence.iterations.inner is the per-outer-iteration cap. interfaceMeasure is the reference integration measure used by the fitting energy. Energy and maximum error answer different questions; neither should silently replace the other as a stopping certificate.

The quality budget is available to spend while fitting. A larger final distortion is not itself failure if it remains below the configured limit. qualityBudgetSatisfied records the sampled checks, not an everywhere-valid certificate for a curved element. Inspect the reason and quality fields even when the target was not reached.

fitting.setMonitor([](const Adaptation::SWIFT::Report& report) {
  // Record or display report.energy, report.geometricSup and quality fields.
});

The monitor receives one snapshot per accepted step and one final report, including early returns. The final snapshot may repeat the last accepted step. It must not modify the problem or throw. getReport() holds the last completed solve. trace additionally enables the examples' diagnostic logs.

Adding metric terms

getMetric() exposes a native additional bilinear form, not the internal scratch problem or a replacement for the canonical fitting and distribution. Both local and global bilinear integrators use Rodin's normal form interface.

fitting.getMetric() += Integral(Dot(u, v));
fitting.getMetric() -= Integral(Dot(u, v));

Added terms are assembled again at each outer iteration. They change the quadratic direction model, not the robust energy or its fitting force. Maintain symmetry, nonnegative curvature and compatible unresolved modes. Arbitrary linear loads are deliberately not exposed: a load without a corresponding fitting energy would invalidate Armijo's descent model.

Implementation layout

All method-specific code lives in Adaptation/SWIFT/: Adapt.h owns the P1 mesh workflow; Problem.h owns outer orchestration; FittingTensor.h and FittingForce.h assemble fitting; Distribution.h assembles distribution; Hinge.h, HingeMetric.h, HingeForce.h and HingeProblem.h define the inner model; LinearSolver.h adapts native linear solvers. Parameters.h and Report.h are the configuration and diagnostic surfaces. Geometry evaluation uses the shared adaptation deformation infrastructure.