Validation

LessCAD runs fast because we rebuilt the geometry kernel and the solvers from the ground up, but fast is only useful if the numbers are right. This page is the verification record: how each solver is tested, and how it performs against known references. For how the engine works, see the Product overview.

Results at a glance

We aim for screening-grade accuracy, within about 15% of standard references, fast enough to catch problems early rather than to replace a certification study. Where a solver has a known limit, we say so plainly.

Overview

Our engine is different in every way. The most important aspects include:

• Implicit Geometry
• Discretization
• Iterative Solvers
• Solution Verification

Implicit Geometry

We chose implicit geometry as our primary representation. This decision shapes every downstream activity, including analysis - and it's the reason our solver stack looks different from conventional FEA.

Two paradigms exist for representing solid geometry in CAD.

We use SDF as our primary geometry kernel. We maintain a B-rep shadow for export to industry-standard formats (STEP, IGES) when the design is ready for fabrication; the SDF carries the authority during authoring and analysis.

Discretization

To solve a partial differential equation on a continuous geometry, we have to discretize - convert the geometry into a finite set of numerical unknowns the solver can iterate on. The discretization choice has cascading consequences for both robustness and speed. Two approaches dominate.

The FCM discrete weak form we solve is:

where:

• $\Omega$: background domain (the FCM grid)
• $\Gamma_N$: Neumann boundary (surfaces where tractions are applied)
• $\alpha(\mathbf{x})$: indicator function (1 inside the part, near zero outside)
• $\mathbf{C}$: elasticity tensor (material properties)
• $u$: displacement field (the unknown we solve for)
• $v$: test function
• $\mathbf{f}$: body force per unit volume (gravity, centrifugal)
• $\mathbf{t}$: surface traction on $\Gamma_N$

The indicator $\alpha(\mathbf{x})$ is 1 inside the part and a small constant $\alpha_\text{out} \approx 10^{-6}$ outside. Outside contributions are nonzero but negligibly small - sufficient to keep the system non-singular on the background grid without polluting the inside-of-part solution.

Immersed Boundary Methods

Three structural problems follow from the boundary not aligning to the grid, and we address each:

• Cut-cell conditioning. A cell that is 99% outside the part contributes a near-zero block to the stiffness matrix, hurting linear solver convergence. We mitigate this with compensated arithmetic in the Krylov inner products and grid-aware preconditioning.
• Dirichlet boundary conditions. A clamped face does not align with grid nodes. We impose Dirichlet BCs weakly via the Nitsche method : the constraint is added as a surface integral that penalizes deviation from the prescribed value, with a stability term that scales with element size and material stiffness. The result is consistent (recovers the strong-form BC as the grid refines) and stable (no ill-conditioning blow-up).
• Surface quadrature. Integrals over cut surfaces require subcell quadrature. We use a marching-cubes-style boundary reconstruction at each cut cell, with Gauss points placed adaptively on the recovered surface.

The tradeoff is real but bounded: FCM gives up the geometric exactness of body-fit meshes in exchange for no meshing pipeline. For our SDF-native authoring model, this is the right exchange. For aspect-ratio-extreme features or razor-thin walls, a coarse FCM grid may underresolve; the solution is a finer grid - the same answer as any FEA code facing the same problem.

Iterative Solvers

Every production solver we ship sits on the same linear-algebra foundation. We made that choice deliberately. Two constraints drove it: the linear systems from our discretization are too large for direct solvers on consumer GPUs, and the work must run entirely in WGSL so it executes in the browser without a CPU fallback.

Why iterative

Direct linear solvers (LU, Cholesky, QR factorizations) produce an exact solution in a finite number of arithmetic operations. They are robust and well-understood. They are also infeasible for our problem class: factorization memory grows as $O(n^{4/3})$ for typical 3D sparsity patterns, exceeding GPU memory budgets at a few hundred thousand unknowns. Production CAD parts produce systems with millions of unknowns. Direct solvers do not finish on the GPUs our users have.

We use iterative solvers instead. They build successive approximations $\mathbf{x}_k$ that converge to the true solution as the iteration count $k$ grows. The termination criterion is not an exact solution but a residual tolerance:

where:

• $A$: system matrix (e.g., stiffness)
• $\mathbf{x}_k$: approximate solution at iteration $k$
• $\mathbf{b}$: right-hand side vector (loads)
• $\varepsilon$: convergence tolerance (typically $10^{-6}$ to $10^{-8}$)

$\varepsilon$ is set by engineering judgment - we use $10^{-6}$ to $10^{-8}$ for production stress and thermal solves. Memory cost is $O(n)$; operations per iteration are matrix-vector products plus inner-products, both of which are highly parallel and map cleanly to GPU kernels.

The cost of going iterative is that convergence depends on the condition number of the system. A well-conditioned system converges in few iterations; an ill-conditioned system can require thousands. FCM systems are ill-conditioned by construction, so preconditioning is where the engineering work happens.

Preconditioning

The conditioning of FCM systems is poor by design - cut cells contribute near-singular blocks. Without preconditioning, even simple problems can require tens of thousands of Krylov iterations to converge, which is prohibitive even on a GPU. We addressed this directly.

A preconditioner $M$ is an approximation of $A$ that is cheap to invert. The iteration is then performed on the preconditioned system:

where:

• $A$: original system matrix
• $M$: preconditioner (an approximation of $A$ that is cheap to invert)
• $\mathbf{x}$: solution vector (the unknown)
• $\mathbf{b}$: right-hand side vector

A good preconditioner satisfies $M \approx A$ in some useful sense while remaining inexpensive enough to apply per iteration. The art is in the balance - a perfect preconditioner ($M = A$) eliminates the outer iteration entirely but is itself the full direct solve we are trying to avoid. A trivial preconditioner ($M = I$) costs nothing per iteration but gives no convergence benefit.

Two preconditioners drive our solver stack:

Block Jacobi. We take $M$ as the block-diagonal of $A$, where each block corresponds to one node's degrees of freedom (three for elasticity, one for thermal). The blocks are inverted in parallel - one workgroup per block. This preconditioner is simple, parallel-friendly, and provides modest convergence benefit. We use it as the baseline preconditioner for our stress and thermal solvers on well-conditioned geometry. Block Jacobi breaks down on systems with strong long-range coupling (slender beams, buckling-class problems) because it cannot represent the low-frequency error components that dominate convergence on such systems.

Geometric multigrid (GMG). We build a hierarchy of grids by successive coarsening of the FCM background grid. Errors are smoothed at each level (via block Jacobi), then restricted to the coarser level, where the residual equation is solved approximately, then prolongated back. The V-cycle refers to the descent-then-ascent pattern through the grid hierarchy. GMG is dramatically more effective than block Jacobi on ill-conditioned systems because it smooths errors at every spatial frequency simultaneously - the coarser grids attack the low-frequency components that block Jacobi cannot reach. The cost is implementation complexity: we wrote restriction and prolongation operators, grid hierarchy bookkeeping, smoother tuning, and block-respecting transfer operators that preserve the displacement DOF structure.

The visualization above shows the down-pass of one V-cycle: a residual pattern appears on the finest 16×16 grid, then each coarser level reveals in turn as the algorithm aggregates that residual upward. Low-frequency components survive the restriction and are visible on the 2×2 grid; high-frequency components average out and are dealt with by the smoother on the fine grid. Each level attacks a different scale of error - that's why GMG converges faster than any single-grid method.

Today, GMG is our production preconditioner for the modal and buckling eigensolvers and an option for stress and thermal PCG. Block Jacobi remains the default for stress and thermal on compact geometry where it is sufficient and cheaper.

Krylov methods

Our outer iteration is a Krylov-subspace method. For symmetric positive definite systems (stress, thermal) we use Preconditioned Conjugate Gradient (PCG). For symmetric generalized eigenvalue problems (modal, buckling) we use LOBPCG. Each is discussed in its solver section below.

The shared property: each Krylov iteration consists of one sparse-matrix-vector product (matrix-free in our implementation - contributions are accumulated cell-by-cell on the fly), a small number of vector inner-products, and the preconditioner apply. All three operations are GPU-native and benefit from the parallelism we already pay for.

Solution Verification

Solution verification asks whether our discretization produces a stable answer that converges as the grid is refined. It is independent of validation - a discretization can be self-consistent (converges to a clean limit) while still solving the wrong problem. ASME V&V 10 names this distinction explicitly.

For our FCM discretization with hex8 trilinear basis functions, theory predicts first-order convergence in the $H^1$ norm:

where:

• $u$: exact (true) solution to the PDE
• $u_h$: discretized (numerical) solution at grid spacing $h$
• $\| \cdot \|_{H^1}$: $H^1$ norm (measures both function and gradient error)
• $h$: characteristic grid spacing
• $p$: polynomial order of basis functions
• $C$: constant depending on geometry, material, and boundary regularity

With $p = 1$ for hex8 trilinear elements. The constant $C$ depends on geometry, material, and boundary regularity. We verify this behavior two ways:

Method of Manufactured Solutions (MMS). We pick a target displacement field $u^*(\mathbf{x})$ - typically a polynomial - and substitute it into the strong-form PDE to derive the body forces $f^*$ and tractions $t^*$ that produce $u^*$ as the exact solution. We then run our solver with those fabricated loads and compare the result against $u^*$. Because $u^*$ is constructed to be the solution, the comparison isolates discretization error from analytic uncertainty. MMS produces order-of-accuracy curves under grid refinement: the observed slope of $\log(error)$ vs. $\log(h)$ should match the theoretical $p$.

Patch tests. We construct small grids with a known analytic answer (uniform tension, simple bending, linear thermal gradient). Our solver must reproduce the answer exactly modulo floating-point round-off. Patch tests catch implementation bugs that would not appear in coarse grids (forgotten signs, transposed indices, BC imposition errors). We run them on every commit as part of the per-PR test suite.

Detailed solution-verification results for each production solver live in the verification archive that ships with each source release.

Linear Stress

Solver

The linear system $K \mathbf{u} = \mathbf{f}$ is symmetric positive definite under physically meaningful BCs. We solve it with Preconditioned Conjugate Gradient (PCG), the canonical Krylov method for this class. The preconditioner is either block Jacobi or block V-cycle multigrid; we pick by problem class - block Jacobi for compact geometry, multigrid for slender or ill-conditioned problems.

Our default convergence criterion is a relative residual tolerance of $10^{-6}$, adjustable per study. We use f32 precision throughout, with compensated summation in the Krylov inner products to maintain accuracy on ill-conditioned systems.

Boundary conditions

Three classes of BC are supported:

• Dirichlet (prescribed displacement): imposed weakly via Nitsche . The user selects a face or surface patch; the BC integral is assembled with a penalty term that scales with element size and material stiffness.
• Neumann (prescribed traction): imposed naturally through the weak form's surface integral on $\Gamma_N$.
• Body forces (gravity, centrifugal): imposed naturally through the weak form's volume integral.

Verification status

Code verification is shipped: we have closed-form analytic references on standard geometries (cantilever, simply-supported beam, plates under uniform loading), patch tests, and selective-DOF BC handling. Our validation against the Code_Aster manual is scaffolded with case SSLV07 (parallelepiped under self-weight) and pending GPU testing integration.

Steady State Thermal

The strong form is the steady heat equation:

where:

• $T$: temperature field (the unknown we solve for)
• $\kappa$: thermal conductivity (material property)
• $q$: volumetric heat source

The discrete weak form follows the same FCM pattern as our stress solver: the background grid carries trilinear hex8 basis functions, $\alpha$-gating handles cells outside the part, and Nitsche imposes Dirichlet BCs on the immersed boundary.

Boundary conditions

• Dirichlet (prescribed temperature): Nitsche weak imposition.
• Neumann (prescribed heat flux): natural integral on $\Gamma_N$.
• Robin (convection): natural integral with the convection coefficient and ambient temperature.

Our thermal solver shares the cross-solver infrastructure described in Iterative Solvers above: PCG with block V-cycle, matrix-free assembly, f32 with compensated arithmetic.

Verification status

Code verification ships analytic references for 1D bar, 2D plate, and 3D brick problems under each BC type, and the convective sphere (a genuinely curved immersed boundary) lands at +1.2% on the temperature rise.

Three known limitations, all open:

• Code_Aster TPLV304 (convective fin) does not pass. The manual gives a tip temperature of 20.329 degC with a band of +/- 0.5; we report 22.1 degC — the fin runs hot. We hold this as a pinned test rather than widening the tolerance. Until it passes, treat convective-fin geometry as unvalidated.
• We cannot hold a fixed temperature on a curved surface. Our fixed-temperature condition pins a thin shell of mesh nodes near the surface rather than the surface itself, and that shell cannot be made thinner than one cell. On a flat face the error is small; on a curved one it is not. Code_Aster TPLV06 (a hollow sphere held at fixed temperature on both curved faces) misses by about half the peak. Fixed-temperature conditions on curved geometry are not validated.
• Flat surfaces are under-integrated, which weakens applied heat fluxes. Our cut-cell integration is accurate on curved boundaries (it converges at the expected second order) but loses about 6-7% of the area on flat, grid-aligned faces, and that error does not shrink with refinement. A heat flux is a density, so it inherits the error directly: a two-layer slab under an applied flux reads 16-28% low.

Energy conservation is checked on every run and holds to machine precision — but we want to be precise about what that does and does not mean. It confirms the solver dropped no term and converged; it is a self-consistency identity, not evidence the answer is right. Both of the failing cases above conserve energy perfectly while being wrong. Accuracy claims on this page come only from the external references, never from that check.

Modal

The governing equation is the generalized eigenvalue problem:

where:

• $K$: elastic stiffness matrix
• $M$: mass matrix
• $\boldsymbol{\phi}_i$: $i$-th mode shape (eigenvector)
• $\omega_i$: $i$-th natural angular frequency (eigenvalue gives $\omega_i^2$)

In context, $K$ is the elastic stiffness matrix, $M$ is the mass matrix, $\boldsymbol{\phi}_i$ is the $i$-th mode shape, and $\omega_i$ is the corresponding natural angular frequency. Engineering quantities of interest are typically the lowest $m$ frequencies - the modes most likely to be excited by external forcing.

Lanczos vs. LOBPCG

Two algorithms are widely used for this class of problem. We originally shipped Lanczos and now use LOBPCG as our production default. The switch was substantiated, not stylistic.

Lanczos shift-invert constructs a Krylov subspace from the operator $(K - \sigma M)^{-1} M$ around a shift $\sigma$. Eigenpairs near $\sigma$ are extracted by projecting onto the Krylov basis and solving a small eigenproblem. The method is robust and well-studied. Its cost structure is unfriendly to our setting: each iteration requires an internal linear solve $(K - \sigma M) \mathbf{v} = M \mathbf{w}$, which itself runs a preconditioned Krylov method to convergence. Lanczos becomes an outer iteration of inner iterations, and the inner-solve tolerance interacts non-trivially with the outer convergence.

LOBPCG (Locally Optimal Block Preconditioned Conjugate Gradient) takes a different approach. It maintains a block of $m$ candidate eigenvectors simultaneously and updates them via a preconditioned conjugate-gradient- style step. There is no inner linear solve. The preconditioner is applied directly to the eigenvalue residual; the block structure lets all $m$ target modes converge in parallel.

LOBPCG fits the infrastructure we already built for three reasons:

• Preconditioner-friendly. The block V-cycle multigrid preconditioner we use for stress is reused directly. No new infrastructure is needed to precondition the eigensolve.
• GPU-friendly. Each LOBPCG iteration is a small fixed sequence of matrix-block products and dense block reductions. The work pattern maps cleanly to WGSL compute kernels.
• Block-aware convergence. All $m$ requested modes converge together, rather than serially. For typical engineering studies with $m = 6$ to $12$, this is a several-fold speedup over Lanczos.

Eigenvalue-locked early stop

For slender geometries, the residual-based convergence criterion can stay loose long after the eigenvalues themselves have stabilized. This matters in practice: the user-facing answer is the natural frequency, not the residual. We monitor eigenvalue stability across iterations and stop when each requested eigenvalue has stayed within a per-solver locked threshold (modal: $10^{-4}$; buckling: $5 \times 10^{-2}$) for several consecutive iterations. The residual-tolerance check still runs in parallel as a safeguard. This is what the lambda-locked badge in the viewport reports.

Verification status

Our code verification uses a slender cantilever against Cowper-corrected Bernoulli-Euler beam theory, computed in-repo. The solver lands at -1.0% on mode 1 at $L/h = 12.5$, and -7.0% at $L/h = 20$. To be precise about what that reference is: it is closed-form beam theory, not a NAFEMS benchmark. We do not hold the NAFEMS suite, and we do not claim it. Bernoulli-Euler itself runs 1-2% high on mode 1 because it ignores shear, so the true solver error is smaller than the headline number. Validation against Code_Aster cases (SDLV401, SDLV124, SDLS03, SDLL101) is scaffolded.

Buckling

The governing equation is another generalized eigenvalue problem:

where:

• $K$: elastic stiffness matrix
• $K_\text{geom}$: geometric stiffness matrix (assembled from a pre-stress field)
• $\boldsymbol{\sigma}$: pre-stress field from a prior stress solve
• $\boldsymbol{\phi}$: buckling mode shape (eigenvector)
• $\lambda$: buckling load multiplier (the smallest positive value is the critical load factor)

In context, $K$ is the elastic stiffness, $K_\text{geom}$ is the geometric stiffness assembled from a pre-stress field $\boldsymbol{\sigma}$, and $\lambda$ is the buckling load multiplier. The smallest positive $\lambda$ is the critical load factor at which the structure first buckles. Practical engineering reporting takes the lowest few buckling modes.

Coupled pipeline

Buckling is not a standalone solve. The pre-stress $\boldsymbol{\sigma}$ must be known before $K_\text{geom}$ can be assembled. Every buckling study we run therefore consists of:

1. A linear stress solve under the applied load, producing a displacement field $\mathbf{u}$ and a stress field $\boldsymbol{\sigma}$.
2. Recovery of $\boldsymbol{\sigma}$ at Gauss points across the FCM grid.
3. Assembly of $K_\text{geom}$ from the recovered stress field.
4. Solution of the generalized eigenvalue problem above.

The eigensolver in step 4 is the same LOBPCG we use for modal, with the same block V-cycle preconditioner.

Known limitation

On slender geometries with $L/h > 20$, our current block V-cycle preconditioner is the binding constraint on accuracy. We underestimate the Euler critical load by roughly a factor of four on test cases at $L/h = 20$. The residual tolerance for the outer eigensolve cannot tighten beyond the preconditioner's smoothing fidelity, which on slender elasticity is bending-mode-limited. Our eigenvalue-locked early-stop mechanism surfaces this honestly: the lambda-locked badge reports "engineering trustworthy" only when stability has converged below the per-solver threshold.

The structural fix is a line smoother that resolves the slender-bending mode directly. It's on our active solver roadmap. Until it ships, treat slender buckling studies as bounded estimates, not precise critical loads.

Verification status

Against the closed-form Euler column ($P_{cr} = \pi^2 E I / 4L^2$), the solver lands at -0.4% at $L/h = 12.5$ and -11.6% at $L/h = 20$.

We hold one PUBLISHED buckling reference, Code*Aster SSLL403 (Greenhill: a column buckling under its own weight, $\lambda_c = 170.701$), and we do not pass it. The measured load factor is roughly two orders of magnitude high. The applied self-weight is exactly right, so the defect is in the eigensolver: the eigenproblem is not invariant under a rescaling of the prestress, though $K + \lambda K*{geom} = 0$ is linear in it. Self-weight buckling of very slender members is therefore NOT trustworthy today, and we would rather say so than widen a tolerance until the row turns green. Point-loaded buckling (the numbers above) is unaffected. SSLV139 (circular plate buckling) is scaffolded. These will surface the $L/h > 20$ limit formally with reference numbers.