diff_soma

A native, high-performance differentiable structural form-finding and multi-objective optimization engine written in Rust with C# Grasshopper integration for Rhino 8 & 9.

---

Overview

diff_soma is a native standalone library designed to facilitate interactive, real-time structural design and gradient-based optimization. By implementing both the Force Density Method (FDM) and an Analytical Adjoint State Method natively in Rust, the engine calculates exact mathematical gradients for multi-objective structural layouts in microseconds. It compiles to a lightweight dynamic library, enabling direct, low-latency integration within desktop CAD platforms like Rhino and Grasshopper without requiring external runtime environments or inter-process communication (IPC) overhead.

---

Technical Pillars

• Non-Linear Forward Pass with Dynamic Self-Weight: Integrates structural dead load updating along the Z-axis based on solved member lengths, resolved via a fast-converging fixed-point iteration loop.
• Analytical Adjoint Gradients: Computes exact gradients of the multi-objective loss with respect to both force densities ($q$) and support coordinates ($x_c$) in $O(E)$ time, utilizing explicit mathematical derivations instead of automatic differentiation packages.
• Block-Triangular Adjoint System: Leverages the block upper-triangular structure of the self-weight Jacobian to solve the coupled 3D adjoint state equations sequentially (back-substitution), bypassing the need for assembling a full $3N_f \times 3N_f$ coupled system.
• Matrix-Free Preconditioned Conjugate Gradient (PCG): Employs a native Jacobi-preconditioned CG solver for symmetric positive-definite linear subsystems. It evaluates matrix-vector products directly from edge connectivity arrays, requiring $O(1)$ auxiliary memory.
• Multi-Objective Loss Formulation: Supports simultaneous optimization of coordinate-wise target geometry matching, projection onto curved parametric target surfaces (Spherical Domes and Hyperbolic Paraboloids), total structural load-path (material volume), discrete Laplacian coordinate smoothing, member length equalization, and quad-panel coplanarization.

---

Mathematical Formulation

1. Forward Pass with Self-Weight

In the Force Density Method (FDM), the nodal equilibrium equations are written as:

$$ \mathbf{D}_f \mathbf{x}_f = \mathbf{p}_x - \mathbf{D}_c \mathbf{x}_c $$

$$ \mathbf{D}_f \mathbf{y}_f = \mathbf{p}_y - \mathbf{D}_c \mathbf{y}_c $$

$$ \mathbf{D}_f \mathbf{z}_f = \mathbf{p}_{z} - \mathbf{D}_c \mathbf{z}_c $$

Where $\mathbf{D}_f$ is the system stiffness matrix computed from the edge force densities $\mathbf{q}$. When accounting for structural self-weight, the load vector $\mathbf{p}_z$ is modeled as a function of the solved coordinates through the current member lengths $L_e$:

$$ p_{i, z} = p_{\text{ext}, i, z} - \frac{1}{2} \sum_{e \in \mathcal{E}_i} w_{\text{self}, e} L_e $$

A fixed-point iteration loop resolves this non-linear dependency in $3$ to $5$ iterations to a coordinate tolerance of $10^{-7}$.

2. Block-Triangular Adjoint Coupling

Let $\mathbf{R} = [\mathbf{R}_x, \mathbf{R}_y, \mathbf{R}_z]^T = \mathbf{0}$ be the nodal residuals. Differentiating $\mathbf{R}$ with respect to the state variables yields the coupled Jacobian $\mathbf{J}$. Due to the dynamic self-weight load depending on the solved coordinates, the transpose of the Jacobian, $\mathbf{J}^T$, exhibits a block upper-triangular structure:

$$ \mathbf{J}^T = \begin{bmatrix} \mathbf{A} & \mathbf{0} & \mathbf{K}_{zx}^T \\ \mathbf{0} & \mathbf{A} & \mathbf{K}_{zy}^T \\ \mathbf{0} & \mathbf{0} & \mathbf{A} + \mathbf{K}_{zz}^T \end{bmatrix} $$

Where $\mathbf{A} = \mathbf{D}_f$ is the symmetric FDM stiffness matrix. By exploiting this block upper-triangular structure, diff_soma solves the coupled 3D adjoint state equations sequentially (back-substitution) rather than assembling a full $3N_f \times 3N_f$ coupled system:

1. Solve the non-symmetric system $(\mathbf{A} + \mathbf{K}_{zz}^T) \boldsymbol{\lambda}_z = \mathbf{g}_z$ using highly optimized LU decomposition.
2. Solve the symmetric systems $\mathbf{A} \boldsymbol{\lambda}_x = \mathbf{g}x - \mathbf{K}{zx}^T \boldsymbol{\lambda}_z$ and $\mathbf{A} \boldsymbol{\lambda}_y = \mathbf{g}y - \mathbf{K}{zy}^T \boldsymbol{\lambda}_z$ using the pre-factored Cholesky decomposition of $\mathbf{A}$.

3. Quad Panel Planarization Loss

To enforce the constructional constraint of flat quadrilateral panels, the coplanarization loss $\mathcal{L}_{\text{planar}}$ is defined for each face $f = (a, b, c, d)$ by the squared volume of the parallelepiped formed by the four vertices:

$$ V_f = (\mathbf{d}_f - \mathbf{a}_f) \cdot ((\mathbf{b}_f - \mathbf{a}_f) \times (\mathbf{c}_f - \mathbf{a}_f)) $$

$$ \mathcal{L}_{\text{planar}} = \frac{1}{2} w_{\text{planar}} \sum_{f} V_f^2 $$

The analytical partial derivatives of the volume $V_f$ are computed efficiently in 3D vector space as:

$$ \frac{\partial V_f}{\partial \mathbf{d}_f} = (\mathbf{b}_f - \mathbf{a}_f) \times (\mathbf{c}_f - \mathbf{a}_f) $$

$$ \frac{\partial V_f}{\partial \mathbf{c}_f} = (\mathbf{d}_f - \mathbf{a}_f) \times (\mathbf{b}_f - \mathbf{a}_f) $$

$$ \frac{\partial V_f}{\partial \mathbf{b}_f} = (\mathbf{c}_f - \mathbf{a}_f) \times (\mathbf{d}_f - \mathbf{a}_f) $$

$$ \frac{\partial V_f}{\partial \mathbf{a}_f} = -\left( \frac{\partial V_f}{\partial \mathbf{b}_f} + \frac{\partial V_f}{\partial \mathbf{c}_f} + \frac{\partial V_f}{\partial \mathbf{d}_f} \right) $$

---

Directory Structure

diff_soma/                # Git Repository Root
├── README.md             # Scientific documentation & FDM formulations
├── LICENSE               # MIT Open-source license
├── .gitignore            # Git exclusion settings (excluding heavy binaries)
│
├── diff_soma/            # Rust Backend Library
│   ├── Cargo.toml
│   └── src/ (lib.rs, fdm.rs, adjoint.rs, optimizer.rs)
│
└── diff_soma_gh/         # C# .NET 8.0 Grasshopper Plugin
    ├── diff_soma_gh.csproj
    ├── diff_soma_ghComponent.cs
    ├── diff_soma_ghInfo.cs
    └── sample/
        └── soma_demo.gh  # Pre-built interactive structural test case

---

Quick Start & Installation

To run the interactive Grasshopper parametric form-finding component inside Rhino 8/9, follow these steps:

1. Download the pre-compiled library files diff_soma_gh.gha and diff_soma.dll from the GitHub Releases page.
2. Open Rhino's Grasshopper components directory. You can easily access this by pressing Win + R on Windows, typing %appdata%\Grasshopper\Libraries, and pressing Enter.
3. Paste both diff_soma_gh.gha and diff_soma.dll directly into that folder.
4. Launch Rhino 8 or Rhino 9, and open Grasshopper.
5. In the Grasshopper tabs, navigate to the newly added Soma tab.
6. Open the pre-built structural interactive demo file located at diff_soma_gh/sample/soma_demo.gh and start tweaking the structural weights, self-weight coefficients, and target geometry parameters to witness real-time multi-objective adjoint form-finding.

---

Verification & Unit Tests

The mathematical correctness of the analytical adjoint gradients has been verified against numerical gradients computed via Central Finite Differences. You can run the precision validation test suite using Cargo:

cargo test

Expected output:

running 1 test
test tests::test_adjoint_precision_vs_finite_differences ... ok

test result: ok. 1 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s

---

Limitations and Future Work

While diff_soma provides a fast, lightweight, and mathematically complete solution for real-time CAD-integrated form-finding, several areas remain open for future development:

• Fully Sparse CSR Direct Solver: Integrating native sparse direct solver packages (such as nalgebra-sparse or faer-sparse) to accelerate non-symmetric system solves in larger structures.
• Material Orthotropy: Modeling orthotropic shell and cable properties to extend the engine's physical accuracy to composite structures.
• Pneumatic Membrane & Cable-net Pre-stress Constraints: Extending the forward and adjoint solvers to handle enclosed volume constraints (pneumatics) to enable real-time air-supported structure form-finding.
• Comparative Performance Benchmarking: Conducting systematic comparative benchmarks against established structural optimization libraries in large-scale scenarios to evaluate runtime and convergence speed.

---

👨‍💻 Author

Developed with passion by Moamin

---

⚖️ License

This project is licensed under the MIT License. See the LICENSE file for details.