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

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

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 and 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

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

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

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

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)

Minimal Working Example

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

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)

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 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. 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 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.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}