⚡ electrical engineering

Efficient finite element quadrature on immersed voxel data

This paper presents an efficient voxel quadrature scheme based on a modified moment-fitting approach that computes custom integration rules to absorb voxel information directly into weights, enabling exact stiffness matrix integration for locally refined finite element spaces while decoupling material resolution from the basis and outperforming traditional pre-integration and composed Gauss-Legendre methods.

Original authors: Matti Christmann, Philipp Kopp, Stefan Kollmannsberger

Published 2026-07-01

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Matti Christmann, Philipp Kopp, Stefan Kollmannsberger

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you have a high-resolution 3D scan of a complex object, like a bone or a piece of machinery, made up of millions of tiny, colored blocks called voxels (think of them as 3D pixels). Each block has its own unique material property, like hardness or density.

Now, you want to run a computer simulation to see how this object reacts to stress, like a bridge holding a car or a bone supporting a body. To do this, engineers usually divide the object into a grid of larger "chunks" (finite elements) to solve the math.

The Problem: The "Mosaic" Dilemma
Here is the catch: Your simulation chunks are often much bigger than the tiny voxels from the scan. One simulation chunk might contain hundreds of different voxels, each with a different material.

To calculate the stress accurately, the computer needs to "taste" every single voxel inside that chunk.

The Old Way (Composed Quadrature): Imagine trying to taste a giant stew by taking a spoonful from every single grain of rice, bean, and carrot inside the pot. It's incredibly accurate, but it takes forever. The computer gets bogged down trying to check every single tiny block.
The "Pre-Integration" Way: To speed things up, engineers used to pre-calculate the "flavor" of every possible chunk size and store it in a library. But this is like having a library of recipes only for specific pot sizes. If you change the size of your pot (refine the mesh) or use a different type of cooking method (change the math basis), your library is useless. You can't adapt on the fly.

The Solution: The "Smart Taster" (Voxel Moment Fitting)
The authors of this paper invented a new method called Voxel Moment Fitting. Think of this as a "Smart Taster."

Instead of tasting every single grain of rice in the pot, the Smart Taster looks at the entire pot, figures out exactly how the different ingredients are mixed, and then creates a custom recipe card (a set of special weights) for that specific pot.

Absorbing the Details: The method takes all the tiny, messy details of the voxels inside a chunk and "absorbs" them directly into the numbers (weights) used for the calculation.
The Magic Result: Once the recipe card is made, the computer only needs to take a few strategic spoonfuls (integration points) to get the exact same result as if it had tasted every single grain.
Adaptability: Because the method calculates the recipe card on the fly based on the current chunk, it works perfectly even if you change the size of the chunks or the math rules in the middle of the simulation. It's like a chef who can instantly adjust the recipe whether you are cooking for 2 people or 200, without needing a pre-written library.

What They Found
The researchers tested this on three scenarios:

A Cube with a Hole: They compared their "Smart Taster" against the old methods. They found that if a simulation chunk contains more than about 6 tiny voxels, their method is significantly faster than the old "pre-calculated library" method, and vastly faster than the "taste-every-grain" method.
An L-Shaped Corner: They used their method to zoom in on a sharp corner where stress is high (adaptive refinement). Because the method doesn't rely on a fixed library, it could zoom in and out smoothly, which the old methods couldn't do easily.
A Cracking Plate (3D): They simulated a 3D plate with a hole and a crack spreading through it. This is a massive calculation.
- Using the old "taste-every-grain" method, this simulation would have taken 118 days to run.
- Using their "Smart Taster," it took only 30 hours.

The Bottom Line
This paper presents a way to make computer simulations of complex, scanned objects much faster without losing accuracy. It does this by turning the messy details of millions of tiny blocks into a few smart numbers, allowing the computer to solve the math quickly and adaptively, even for very large and complex 3D problems.

Technical Summary: Efficient Finite Element Quadrature on Immersed Voxel Data

Problem Statement
High-resolution computed tomography (CT) scans provide detailed voxel-wise information regarding geometry and material properties. While direct numerical simulation treating every voxel as a finite element is possible, it often results in an excessive number of unknowns (billions). Coarsening the mesh to reduce degrees of freedom introduces a significant challenge: accurate numerical quadrature. When a single finite element encompasses multiple voxels with heterogeneous material properties (e.g., from CT grey values), the integrand of the stiffness matrix becomes discontinuous. Standard Gauss-Legendre quadrature, which assumes smooth polynomial integrands, yields large integration errors in this context.

The standard remedy, "composed integration," subdivides the element integral at voxel boundaries, effectively performing integration on every voxel. While accurate, this is computationally expensive due to the massive number of integration points. An alternative, "pre-integration," computes voxel-specific sub-stiffness matrices once and assembles them later. However, pre-integration has critical limitations: it cannot be easily adapted when the number of voxels per element varies, when polynomial degrees change, or, most importantly, when adaptive mesh refinement is employed. In adaptive schemes (such as multi-level $hp$-methods), the basis functions vary across refinement levels, rendering pre-computed sub-matrices invalid.

Methodology
The authors propose a voxel-based moment fitting scheme to address these limitations. This method extends the moment fitting approach to immersed voxel geometries. The core idea is to derive a custom integration rule for each finite element that absorbs the discontinuous material distribution directly into the integration weights, rather than the integration points.

Mathematical Formulation: The element stiffness integral is split into a discontinuous material function $\alpha(x)$ and a smooth part $h(x)$ . Instead of integrating the product directly, the method solves a system of linear equations to determine custom weights $w_j$ for a fixed set of integration points $x_i$ .
Lagrange Basis: By choosing Lagrange basis functions defined at the integration points, the system of moment equations diagonalizes. The weights are computed as moments: $w_j = \int_{\Omega_e} \alpha(x) l_j(x) d\Omega$ .
Voxel Integration: The right-hand side of the moment equation is evaluated using composed quadrature on the voxel level. Crucially, the authors introduce a pre-integrated scheme (Algorithm 2) where the geometric part of the moment ( $\int l_j(x) dx$ over a voxel) is pre-computed into a matrix $S$ . The final weights for any element are then obtained by a simple contraction of $S$ with the element's specific material values $\alpha_k$ . This avoids re-evaluating integrals over voxels for every element.
Adaptivity: Because the weights are computed based on the material distribution and the specific basis functions of the current element, the method is independent of the global mesh topology. This allows it to function seamlessly with multi-level $hp$-adaptivity, where basis functions change across refinement levels.

Key Contributions

Basis-Independent Quadrature: The method enables exact integration of stiffness matrices for elements with any polynomial degree and arbitrary voxel layouts without requiring pre-integrated sub-matrices. This unlocks adaptive refinement strategies (like multi-level $hp$) that were previously incompatible with voxel-based pre-integration.
Computational Efficiency: The scheme absorbs per-voxel material information into custom quadrature weights. This eliminates the need for a costly loop over all voxels during the assembly of the global stiffness matrix, a bottleneck in composed quadrature.
Pre-integrated Moment Matrix: The introduction of the pre-computed matrix $S$ significantly reduces the computational cost of determining weights, decoupling the material resolution from the finite element basis.

Results
The paper validates the method through three numerical examples:

Scaling Analysis (Heterogeneous Unit Cube): A comparison of composed quadrature, pre-integration, and voxel moment fitting on a cube with $768^3$ $76 8^{3}$ voxels.
- Accuracy: Voxel moment fitting yields identical strain energy to composed quadrature, confirming mathematical equivalence.
- Performance: For elements containing more than $\approx 6$ voxels per side length (at $p=1$ ), voxel moment fitting outperforms pre-integration. At high voxel counts (e.g., 384 voxels per side), it achieves a runtime reduction of 99.6% compared to composed quadrature and 89.7% compared to pre-integration.
- Break-even Point: The efficiency crossover depends on the polynomial degree $p$ ; higher $p$ shifts the break-even point toward elements with more voxels.
Adaptive L-shaped Domain (2D): The method is applied to a plane-stress problem with a re-entrant corner singularity using multi-level $hp$-adaptivity.
- The scheme successfully resolves the singularity and the geometric stair-stepping of the voxelized boundary.
- The study quantifies the interaction between discretization error and geometric error, showing that adaptive refinement should stop once elements are refined below the voxel size, as further refinement only resolves the voxelized geometry rather than the physical domain.
Phase-Field Fracture (3D): A simulation of a notched plate with a hole and spatially varying Young's modulus.
- The method handles the coupling of mechanical equilibrium and phase-field evolution. Although the degradation function $g(s)$ raises the integrand degree (leading to under-integration of the stiffness), the results show that limit loads and crack paths are largely unaffected.
- Feasibility: A full simulation with 50 load steps took 30 hours using voxel moment fitting. The authors estimate that using composed quadrature would have required 118 days, rendering the problem computationally infeasible without this scheme.

Significance
The paper claims that voxel-based moment fitting provides a robust and efficient pathway for high-fidelity simulations on CT-derived data. Its primary significance lies in decoupling material resolution from the finite element basis, thereby enabling adaptive refinement strategies that were previously impossible with pre-integration techniques. By reducing the computational cost of quadrature by orders of magnitude compared to composed integration, the method makes large-scale, 3D, heterogeneous fracture simulations feasible on standard hardware. The authors emphasize that the method requires no changes to standard finite element codes beyond the weight calculation, making it a practical addition to existing solvers.

Technical Summary: Efficient Finite Element Quadrature on Immersed Voxel Data

More like this