A ROM-based BDDC solver for unfitted p-FEM level-set-based lattice structures

This paper presents a scalable ROM-enhanced BDDC solver that enables fast, high-accuracy simulations of complex, graded lattice structures using unfitted p-FEM without relying on homogenization or multiscale assumptions.

The Gist

The Big Picture: Building a Digital Lego Castle

Imagine you are an engineer trying to design a super-lightweight, incredibly strong bridge or a custom wrench. Instead of using solid metal, you want to build it out of a complex, honeycomb-like lattice (like a bird's bone or a sponge).

The problem? These structures are made of thousands of tiny, unique cells. If you try to simulate how they bend or break on a computer, the math gets so heavy that it would take a supercomputer days to solve it. Traditional methods often try to "cheat" by pretending the whole thing is a uniform block of material, but that fails when the design gets fancy and changes from one spot to another.

This paper presents a new "super-fast" way to simulate these complex structures on a regular laptop, solving a problem with 17,000 tiny cells in just 30 seconds.

Here is how they did it, broken down into three main tricks:

---

1. The "One-Size-Fits-All" Grid (Unfitted p-FEM)

The Problem: Usually, to simulate a shape, you have to draw a mesh (a grid of tiny triangles) that fits perfectly inside the shape. If the shape is a weird, cut-out circle, you have to cut the triangles to fit the curve. If you have 17,000 different shapes, you have to cut 17,000 different sets of triangles. That's a nightmare.

The Solution: Imagine you have a giant, transparent sheet of graph paper (the grid). You place your weird shapes on top of it. You don't cut the paper; you just tell the computer, "Only count the parts of the paper that are inside the shape."

• The Analogy: Think of it like a cookie cutter on a sheet of dough. You don't need to reshape the dough; you just ignore the parts outside the cutter.
• The "High-Order" Twist: Instead of using tiny, simple triangles (like low-resolution pixels), they use giant, super-smooth, high-definition "pixels" (polynomials). This means they can get a very accurate picture of the curve using very few pieces, keeping the computer's workload low.

2. The "Divide and Conquer" Strategy (BDDC)

The Problem: Even with the grid trick, solving the math for 17,000 cells all at once is like trying to solve a 17,000-piece jigsaw puzzle by looking at every single piece at the same time.

The Solution: They split the puzzle into 17,000 tiny, separate puzzles (one for each cell).

• The Analogy: Imagine a massive library. Instead of one librarian trying to find a book in the whole building, they give each section of the library its own librarian.
• The "Coarse" Step: The local librarians solve their small sections quickly. But then, they need to make sure the books match up at the borders. A "Head Librarian" (the coarse solver) steps in to balance everything out and ensure the whole library makes sense. This method (BDDC) is famous for being able to handle huge problems efficiently by keeping the local work independent.

3. The "Cheat Sheet" (Reduced Order Modeling & ROM)

The Problem: This is the biggest bottleneck. Even with the grid and the divide-and-conquer tricks, the computer still has to do a massive amount of heavy math (integration) for every single cell to figure out how stiff it is. Doing this 17,000 times is slow.

The Solution: They realized that while every cell is in a different spot, they are all made of the same "recipe" (the same level-set function, like a Schwarz Diamond).

• The Analogy: Imagine you are baking 17,000 cookies. You don't need to weigh the flour and sugar for every single cookie from scratch. You bake a few "training" cookies, measure exactly how much flour and sugar you used, and write down a "Cheat Sheet" (a Reduced Order Model).
• How it works:
  1. Offline Training: They do the hard math once for a few sample shapes and store the results in a database.
  2. Online Speed: When they need to simulate the 17,000 cells, they don't do the heavy math. They just look up the "Cheat Sheet" and do a quick multiplication.
  3. The Magic: This "Cheat Sheet" is built using a technique called MDEIM (Matrix Discrete Empirical Interpolation Method). It's like a smart compression algorithm that remembers the most important parts of the math so it doesn't have to recalculate the boring stuff.

The "Stabilizer" (The Safety Net)

Because they are using these "Cheat Sheets" and ignoring parts of the grid, the math can sometimes get a little wobbly or unstable.

• The Fix: They added a tiny "stabilization term." Think of this like adding a little bit of glue to the edges of your puzzle pieces. It doesn't change the picture, but it stops the pieces from wobbling around, ensuring the computer doesn't crash or give a wrong answer.

The Result: A Laptop Supercomputer

By combining these three ideas:

1. Smart Grids (ignoring the messy edges),
2. Divide and Conquer (splitting the work), and
3. Cheat Sheets (skipping the heavy math),

The authors managed to simulate a complex 2D lattice structure with 17,000 unique cells on a standard laptop in 30 seconds.

Why does this matter?
In the past, engineers had to simplify their designs or wait days for a supercomputer to check if a new lightweight material would work. Now, they can design complex, graded lattices (where the material gets stronger or weaker in specific spots) and test them instantly. This opens the door to creating ultra-lightweight, super-strong materials for aerospace, medical implants, and sports equipment much faster than ever before.

Technical Summary

1. Problem Statement

The paper addresses the computational challenge of simulating large-scale lattice structures (architected materials) with complex, graded geometries.

• Geometric Complexity: These structures are often defined by level-set functions (e.g., Triply Periodic Minimal Surfaces or TPMS) rather than explicit parametric mappings. This results in "unfitted" meshes where the geometry cuts through a background grid, creating cut elements.
• Limitations of Existing Methods:
  • Homogenization/Multiscale methods: Rely on assumptions of scale separation and periodicity, which fail for graded lattices or slender structures with few cells in the thickness direction.
  • Standard FEM: Direct full-scale 3D simulations are computationally prohibitive due to the massive number of degrees of freedom (DoFs) and the complexity of integrating over cut elements.
  • Existing Accelerated Solvers: Many fast solvers rely on conforming discretizations and parametric morphing to exploit cell similarities via Reduced Order Modeling (ROM). These cannot handle the arbitrary topologies and non-parameterizable geometries typical of level-set-based cellular structures.

2. Methodology

The authors propose a novel solver combining Unfitted High-Order Finite Element Methods (p-FEM), Balanced Domain Decomposition by Constraints (BDDC), and Reduced Order Modeling (ROM).

A. Geometric Modeling & Discretization

• Level-Set Definition: Cells are defined by mapping a parametric unit square to physical space and trimming it using a level-set function $\phi < \mu$. The threshold $\mu$ can vary linearly across the cell, allowing for graded designs.
• Unfitted p-FEM: The domain is discretized using a single high-order element per cell (e.g., degree $p=8$) defined on the untrimmed parametric domain.
  • Advantage: Basis functions remain active regardless of the trimming, ensuring the stiffness matrix structure is continuous with respect to geometric changes. This is crucial for ROM.
  • Stabilization: An $\alpha$-stabilization term is added to the weak form to handle the ill-conditioning of cut elements, interpreted as immersing the cut region in a soft material.

B. Domain Decomposition (BDDC)

• Subdomain Definition: Each lattice cell corresponds to one subdomain.
• Solver Strategy: The global system is solved using the Preconditioned Conjugate Gradient (PCG) method with a BDDC preconditioner.
  • Coarse Space: Includes cross-point values, edge averages, and edge moments to ensure scalability.
  • Parallelism: Local Neumann problems are solved independently (parallel), while a global coarse problem balances the corrections.

C. Acceleration via ROM and Fast Assembly

The primary bottleneck is the assembly of stiffness matrices for cut elements, which requires expensive quadrature. The authors accelerate this via:

1. Fast Assembly Technique: The geometric mapping contributions are separated from the trimmed domain integration. The integrands are approximated by polynomials, allowing the stiffness matrix to be computed via tensor contractions rather than full numerical integration.
2. Matrix Discrete Empirical Interpolation Method (MDEIM):
   • A Reduced Order Model (ROM) is built offline to approximate the integrals of the fast assembly tensors.
   • Clustering: The parameter space (threshold values) is divided into clusters. A specific ROM is trained for each cluster using Latin Hypercube Sampling.
   • Online Phase: During simulation, the ROM surrogate replaces the expensive quadrature for cut elements. The ROM is valid for any geometric mapping, provided the level-set function remains the same.
3. Stabilization for ROM: A specific stabilization term is introduced to ensure the BDDC solver remains scalable when using the ROM approximation, preventing iteration counts from exploding.

3. Key Contributions

• Novel Solver Framework: Integration of BDDC with unfitted p-FEM and ROM, specifically tailored for non-periodic, level-set-defined lattice structures.
• ROM for Unfitted Methods: Development of a ROM strategy that works with unfitted discretizations by exploiting the continuity of the basis functions across trimming variations, bypassing the need for conforming meshes.
• Scalability with Graded Designs: The method avoids homogenization assumptions, allowing for the simulation of structures with varying cell geometries and topologies (graded lattices).
• Efficiency: The combination of fast assembly and ROM reduces the assembly time significantly, making large-scale simulations feasible on standard hardware.

4. Numerical Results & Performance

The method was validated on a MacBook Pro M4 Pro (14 cores, 48GB RAM) using 8 parallel cores.

• Accuracy:
  • High-order p-FEM ($p=8$) achieves accuracy comparable to refined CutFEM meshes but with significantly fewer DoFs.
  • The ROM error is controlled and balanced with the fast assembly interpolation error.
• Scalability:
  • The number of solver iterations remains asymptotically bounded (converging to ~17 iterations) as the number of subdomains increases, provided the ratio of subdomain size to mesh size ($H/h$) is constant.
• Speed:
  • 17,000 Cells: A complex 2D problem with over 17,000 cells of varying geometry was solved in approximately 30 seconds.
  • Comparison: The proposed method is ~10x faster than Cholesky decomposition and ~100x faster than SOR-preconditioned Conjugate Gradient for these problems.
• Applications:
  • Sandwich Wing: A 4,000-cell wing profile with a porous core and solid skin, solving under aerodynamic loads in 18 iterations.
  • Lattice Wrench: A non-tensor-product layout with 90 subdomains, demonstrating robustness for irregular geometries.

5. Significance

This work represents a significant leap in the simulation of architected materials. By removing the reliance on homogenization and periodicity, it enables the design and analysis of graded, complex lattice structures that are increasingly common in additive manufacturing.

• Design Optimization: The speed of the solver (seconds for thousands of cells) makes it a viable candidate for integration into topology optimization loops, where thousands of forward solves are required.
• Gradient Computation: The ROM surrogate facilitates the efficient computation of gradients with respect to design parameters.
• Future Work: The authors plan to extend this to 3D problems (addressing the "curse of dimensionality" in ROM), hyperelastic materials, and plate simulations.

In summary, the paper presents a highly efficient, scalable, and accurate computational framework that bridges the gap between complex geometric descriptions (level-sets) and large-scale structural analysis, overcoming the limitations of traditional multiscale and full-scale FEM approaches.