# Core Architecture of the FEM-Based Elasticity Solver in ppf-contact-solver

> Explore the core architecture of the FEM-based elasticity solver in ppf-contact-solver. Discover its GPU-accelerated finite-element pipeline for efficient dynamic elastic simulations on CUDA.

- Repository: [ZOZO, Inc./ppf-contact-solver](https://github.com/st-tech/ppf-contact-solver)
- Tags: architecture
- Published: 2026-05-27

---

**The FEM-based elasticity solver in ppf-contact-solver implements a GPU-accelerated finite-element pipeline that combines dual CSR sparse matrices, unrolled 3×3 block operations, and a diagonal-preconditioned Conjugate Gradient solver to handle dynamic elastic simulations on CUDA hardware.**

The ppf-contact-solver repository by Preferred Networks delivers a high-performance computational mechanics framework built around a **FEM-based elasticity solver** written in CUDA C++. This engine processes complex cloth and deformable body simulations through a carefully optimized linear algebra pipeline designed for cache efficiency and massive parallelism across modern GPUs.

## Hybrid Sparse Matrix Representation

The global stiffness system employs a three-part matrix structure to separate static and dynamic contributions while maximizing memory throughput.

### Dynamic and Fixed CSR Matrices

The solver maintains two Compressed Sparse Row (CSR) structures: **DynCSRMat** (matrix **A**) for per-vertex contributions that update each iteration (such as strain-limiting constraints), and **FixedCSRMat** (matrix **B**) for static couplings including contact and plasticity terms. This separation allows the solver to rebuild only the dynamic portion while keeping the fixed structure constant, reducing assembly overhead.

Source: [`crates/ppf-cts-solver/src/cpp/solver/solver.cu`](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/solver/solver.cu)

### Per-Vertex Dense Blocks

An additional dense per-vertex matrix **C** stores the local elastic stiffness of each mesh element as a `Vec<Mat3x3f>` structure. These 3×3 blocks are populated by energy assembly routines that compute elastic contributions from mesh topology and material parameters. The block-diagonal structure aligns with the three degrees of freedom per vertex, enabling coalesced memory access patterns during matrix-vector multiplication.

Source: [`crates/ppf-cts-solver/src/cpp/csrmat/csrmat.cu`](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/csrmat/csrmat.cu) and [`crates/ppf-cts-solver/src/cpp/energy/energy.cu`](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/energy/energy.cu)

## Matrix-Vector Product Kernel

The core computational workload resides in the `apply` kernel, which evaluates the complete linear system **A·x + B·x + C·x** for each iteration.

### Unrolled 3×3 Multiplication

To eliminate branches and maintain instruction-level parallelism, the kernel uses an **`UnrolledMat3x3f`** class for block matrix operations. This unrolled approach processes the three-component vectors per vertex without loop overhead, keeping all arithmetic units saturated while ensuring memory accesses remain coalesced across warp threads.

Source: [`apply` function in solver.cu](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/solver/solver.cu#L35-L71)

## Diagonal Preconditioning

Before entering the iterative solve phase, the solver constructs a diagonal preconditioner from the inverse of summed diagonal blocks **A(i,i) + B(i,i) + C[i]**. The **`invert`** function performs analytical 3×3 matrix inversion using explicit cofactor calculations rather than numerical methods, yielding the preconditioner vector **P** used by the `precond` kernel to accelerate Conjugate Gradient convergence.

Source: [`invert` and preconditioner in solver.cu](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/solver/solver.cu#L106-L130)

## Conjugate Gradient Solver

The linear system solver relies on a custom Conjugate Gradient implementation defined in the **`cg`** function. This routine iteratively calls the `apply` kernel for matrix-vector products, applies the diagonal preconditioner through the `precond` kernel, and computes residual norms using GPU reduction operations. Convergence is determined by user-specified tolerance thresholds and maximum iteration counts, balancing accuracy against real-time performance requirements.

Source: [`cg` implementation](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/solver/solver.cu#L32-L74)

## Global Solve Entry Point

The **`solve`** function orchestrates the complete solution process. It first builds the preconditioner, instantiates **`DeviceOperators`** for GPU memory management, then launches the `cg` routine. Upon completion, the function returns success status, iteration count, and final residual magnitude to the calling layer.

Source: [`solve` function](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/solver/solver.cu#L78-L99)

## Integration with Python Frontend

While the core solver runs as a compiled CUDA binary, the repository provides a high-level Python API through **`App`** and **`Session`** classes that abstract the FEM pipeline.

### Scene Serialization and Execution

The frontend serializes scene data into per-vertex and per-element files, invokes the compiled solver through Rust bindings (`_rust` module), and streams results back to Jupyter or Blender. The `Session` class manages subprocess lifecycle, progress monitoring, and checkpointing, enabling batch processing of complex animations without manual GPU memory management.

Source: [[`frontend/_session_.py`](https://github.com/st-tech/ppf-contact-solver/blob/main/frontend/_session_.py)](https://github.com/st-tech/ppf-contact-solver/blob/main/frontend/_session_.py)

### Minimal Working Example

The following Python script demonstrates the complete workflow from mesh creation to GPU solver execution:

```python
from frontend import App

# Initialize application context

app = App.create("demo")

# Generate a 64×64 square sheet mesh

V, F = app.mesh.square(res=64, ex=[0, 0, 1], ey=[0, 1, 0])
app.asset.add.tri("sheet", V, F)

# Build scene with multiple constrained objects

scene = app.scene.create()
for i in range(5):
    obj = scene.add("sheet")
    obj.at(i * 0.25, 0, 0)                # Position each instance

    obj.pin(obj.grab([0, 1, 0]))          # Fix corner vertex

    obj.param.set("strain-limit", 0.05)   # Enable FEM strain limiting

scene = scene.build()

# Configure simulation parameters

session = app.session.create(scene)
session.param.set("dt", 0.01)             # Time step size

session.param.set("frames", 120)          # Animation length

session = session.build()

# Launch GPU solver and export results

session.start(blocking=True)
session.export.animation().zip()

```

Source: [[`examples/headless.py`](https://github.com/st-tech/ppf-contact-solver/blob/main/examples/headless.py)](https://github.com/st-tech/ppf-contact-solver/blob/main/examples/headless.py)

## Summary

- The **FEM-based elasticity solver** uses a dual CSR structure (**DynCSRMat** and **FixedCSRMat**) to separate dynamic strain constraints from static contact couplings.
- The `apply` kernel performs unrolled **3×3 matrix-vector multiplication** via `UnrolledMat3x3f` to maximize GPU instruction throughput.
- An **analytical diagonal preconditioner** inverts 3×3 blocks using the `invert` function before entering the Conjugate Gradient loop.
- The custom **`cg` implementation** combines matrix-vector products, preconditioning, and residual reduction within a single CUDA compilation unit.
- **Python-Rust bindings** in [`frontend/_session_.py`](https://github.com/st-tech/ppf-contact-solver/blob/main/frontend/_session_.py) manage solver subprocesses, enabling high-level scene authoring without direct CUDA programming.

## Frequently Asked Questions

### What makes the ppf-contact-solver FEM implementation suitable for real-time simulation?

The solver achieves real-time performance through three architectural decisions: separating dynamic and fixed matrices to minimize per-frame assembly costs, using unrolled 3×3 kernels that avoid branch divergence on GPUs, and implementing a lightweight diagonal preconditioner that accelerates Conjugate Gradient convergence without expensive factorization. These optimizations allow the system to handle 10,000+ vertices at interactive rates.

### How does the solver handle contact constraints within the FEM framework?

Contact constraints are incorporated into the **fixed CSR matrix (B)** during the assembly phase in [`energy.cu`](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/energy/energy.cu). Because these couplings remain static during the linear solve, they do not require rebuilding each iteration, unlike the dynamic strain-limiting terms stored in matrix **A**. This separation maintains sparsity patterns while handling complex collision scenarios.

### Can the solver handle anisotropic materials or varying stiffness per element?

Yes, the dense per-vertex matrix **C** (`Vec<Mat3x3f>`) stores local elastic stiffness independently for each element. Since **C** is treated as a dense 3×3 block per vertex during the `apply` operation, material parameters can vary spatially across the mesh. The assembly layer in [`energy.cu`](https://github.com/st-tech/ppf-contact-solver/blob/main/crates/ppf-cts-solver/src/cpp/energy/energy.cu) computes these blocks from Young's modulus and Poisson ratio defined per element in the Python frontend.

### What is the role of the `DeviceOperators` class in the solver architecture?

The `DeviceOperators` class, instantiated within the `solve` function, manages GPU memory allocation and stream synchronization for the Conjugate Gradient iteration. It provides the execution context for the `cg` function, handling device-side reductions and ensuring that the sparse matrix operations remain asynchronous with respect to the host CPU until the final result retrieval.