Factorized FEM, per-point quadrature data, and variable-length summations (1.1.1)

.github/workflows/ci.yml

@@ -11,11 +11,11 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-24.04, macos-14, windows-2022]
+        os: [ubuntu-26.04, macos-26, windows-2025]
 
     steps:
       - name: Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v7
 
       - name: Install libomp (macOS)
         if: runner.os == 'macOS'

.gitignore

@@ -1,3 +1,4 @@
 __pycache__
 build/
 build/docs
+docs/gallery/

changelog.md

@@ -1,5 +1,21 @@
 # Changelog
 
+## 1.1.1
+
+- Added factorized FEM kernels with exact derivatives for node-local position maps,
+  including rectangular deformation gradients for lines and membranes.
+- Added per-point quadrature data and variable-length summations per element.
+- Added symbolic substitution for scalars, vectors, and matrices.
+- Added per-potential energy and gradient breakdowns.
+- Direct LLT solves now build the sparse matrix directly from the blocked Hessian
+  and reuse its storage across solves.
+- Memoized expression evaluation over shared subexpressions.
+- `SYMX_ENABLE_MKL` now defaults to `AUTO` and uses oneMKL/PARDISO when found.
+- Added `SYMX_USE_SYSTEM_FMT` (default `OFF`); the bundled fmt is used unless enabled.
+- Updated bundled fmt to 12.2.0.
+- Added the "Advanced FEM and Integration" docs page and factorized FEM options
+  in the rubber extrusion example.
+
 ## 1.1.0
 
 - Added stencil-based differentiable simulation with adjoint gradients,

docs/source/_ext/symx_pygments.py

@@ -64,7 +64,8 @@ def __init__(self, **options):
                     "GlobalPotential", "Verbosity", "NewtonsMethod", "ProjectionToPD",
                     "SolverReturn", "LabelledConnectivity", "numeric_limits", "LinearSolver",
                     "Workspace", "Compiled", "DiffCache", "Compilation", "BlockedSparseMatrix", "bsm",
-                    "Ordering", "RowMajor", "Tet4", "Hex8", "Tet10", "Hex27", "Element", "VectorXd", "symx"
+                    "Ordering", "RowMajor", "Tet4", "Hex8", "Tet10", "Hex27", "Element", "VectorXd", "symx",
+                    "FemEnergySplit"
                 ],
                 funcs=[
                     "PLACEHOLDER",

docs/source/advanced_integration.md

@@ -0,0 +1,208 @@
+# Advanced FEM and Integration
+
+The [FEM integration](fem_integration.md) and [symbol-data maps](symbol_data_maps.md) pages introduced scalar potential definitions and shared, fixed integration rules.
+This page covers:
+
+- Factorized FEM potentials for efficient Hessian generation
+- Element-specific data at integration points
+- Variable summation lengths per element
+
+
+## Factorized FEM Potentials
+
+For high-order elements, differentiating the complete scalar energy with respect to every nodal coordinate can produce large symbolic expressions and generated kernels.
+Declaring the material density separately from the deformation gradient lets SymX differentiate the smaller material expression and apply the remaining matrix products in compact loops.
+This greatly reduces symbolic differentiation, code generation, and compilation time for high-order elements, with a modest gain in evaluation speed.
+The element energy, gradient, and Hessian are the same as those of the standard definition, up to round-off.
+
+At one integration point, collect the nodal positions in $\mathbf x\in\mathbb R^{n\times3}$.
+The deformation gradient is linear in these positions:
+
+$$
+\mathbf F = \mathbf x^\mathsf T\mathbf B,
+\qquad E_q = \psi(\mathbf F),
+\qquad \psi(\mathbf F) = w_q\det(\mathbf J_{\mathbf X})\,\Psi(\mathbf F),
+$$
+
+where $\mathbf B\in\mathbb R^{n\times r}$ is independent of the current positions.
+For standard isoparametric elements, $\mathbf B$ contains the shape function gradients in rest coordinates, and $r$ is the reference dimension: 3 for solids, 2 for membranes and shells, and 1 for lines.
+For $r<3$, $\mathbf F\in\mathbb R^{3\times r}$ is rectangular, and the material density is written in terms of it, for example through $\mathbf C=\mathbf F^\mathsf T\mathbf F$, and $\det(\mathbf J_{\mathbf X})$ is replaced by the corresponding area or length element.
+Writing the coordinates and deformation gradient as vectors $\mathbf q$ and $\mathbf f$, this is a linear map $\mathbf f=\mathbf L\mathbf q$.
+The chain rule gives
+
+$$
+\nabla_{\mathbf q}E_q = \mathbf L^\mathsf T\nabla_{\mathbf f}\psi,
+\qquad
+\nabla^2_{\mathbf q}E_q = \mathbf L^\mathsf T\bigl(\nabla^2_{\mathbf f}\psi\bigr)\mathbf L.
+$$
+
+For Newton's method, SymX generates the material gradient and Hessian symbolically, then the factorized kernel forms the full nodal gradient and Hessian from them.
+The material density $\Psi$ can be nonlinear.
+
+Consider the Stable Neo-Hookean definition from [Newton's method](newtons_method.md).
+Here `data.element_type` has type `FEM_Element`, and the flat connectivity has `data.stride` entries per element.
+The standard scalar definition is:
+
+```cpp
+G->add_potential("stable_neo_hookean_" + get_name(data.element_type),
+    data.connectivity, data.stride,
+    [&](MappedWorkspace<double>& mws, Element& elem)
+    {
+        // Create symbols from data
+        std::vector<Vector> x = mws.make_vectors(data.x, elem);
+        std::vector<Vector> X = mws.make_vectors(data.X, elem);
+        Scalar E = mws.make_scalar(data.youngs_modulus);
+        Scalar nu = mws.make_scalar(data.poissons_ratio);
+
+        return fem_integrator(mws, data.element_type,
+            [&](Scalar& w, Vector& xi)
+            {
+                Matrix Jx = fem_jacobian(data.element_type, x, xi);
+                Matrix JX = fem_jacobian(data.element_type, X, xi);
+                Matrix F = Jx * JX.inv();
+                return stable_neohookean_energy_density(F, E, nu) * JX.det() * w;
+            }
+        );
+    }
+);
+```
+
+With `add_potential_fem()`, return a `FemEnergySplit` containing the nodal position matrix, the coefficients `B`, and the weighted material density:
+
+```cpp
+G->add_potential_fem("stable_neo_hookean_" + get_name(data.element_type),
+    data.connectivity, data.stride,
+    [&](MappedWorkspace<double>& mws, Element& elem)
+    {
+        std::vector<Vector> x = mws.make_vectors(data.x, elem);
+        std::vector<Vector> X = mws.make_vectors(data.X, elem);
+        Scalar E = mws.make_scalar(data.youngs_modulus);
+        Scalar nu = mws.make_scalar(data.poissons_ratio);
+        Matrix x_m(collect_scalars(x), {int(x.size()), 3});
+
+        auto rule = get_integration_rule(data.element_type);
+        return mws.add_for_each(rule, [&](Vector& ip)
+        {
+            Scalar w = ip[0];
+            Vector xi({ip[1], ip[2], ip[3]});
+            Matrix JX = fem_jacobian(data.element_type, X, xi);
+            Matrix B = fem_deformation_gradient_B(data.element_type, xi, JX);
+            auto psi = [=](const Matrix& F) {
+                return stable_neohookean_energy_density(F, E, nu) * JX.det() * w;
+            };
+            return FemEnergySplit{ x_m, B, psi };
+        });
+    }
+);
+```
+
+Both versions use the same integration rule and define the same discrete energy.
+Capture local symbols by value in `psi`, since SymX calls it after the definition callback returns.
+Its dependence on the current positions must be through `F`; `B` must be independent of the state DoFs.
+As with `add_potential()`, the callback can also return a `std::pair<FemEnergySplit, Scalar>` to add a potential-level condition.
+The material density itself may contain branches.
+
+In either version, register the current positions as DoFs:
+
+```cpp
+G->add_dof(data.x);
+```
+
+The [rubber extrusion showcase](rubber-extrusion-example) includes standard, indexed-quadrature, and factorized definitions in [`rubber_extrusion.cpp`](../../examples/rubber_extrusion.cpp), using the material density from [elasticity_potentials.h](../../examples/elasticity_potentials.h).
+
+Factorized SIMD kernels are not provided. The full symbolic energy remains available to the first-order optimizer and adjoint solver.
+
+
+### Mapped Node Positions
+
+The previous example uses nodal positions directly as DoFs.
+Each node's position can also be a nonlinear expression in that node's three DoFs.
+SymX detects this mapping in `FemEnergySplit::x` and handles its derivatives automatically; `B` remains independent of the state DoFs.
+
+Collecting these mapped positions into $\mathbf q=\boldsymbol\chi(\mathbf u)$, let $\mathbf J_\chi=\partial\mathbf q/\partial\mathbf u$.
+Writing the position gradient and Hessian as $\mathbf g_q=\nabla_{\mathbf q}E_q$ and $\mathbf H_q=\nabla^2_{\mathbf q}E_q$, the derivatives with respect to the optimization variables are
+
+$$
+\nabla_{\mathbf u}E_q = \mathbf J_\chi^\mathsf T\mathbf g_q,
+\qquad
+\nabla^2_{\mathbf u}E_q = \mathbf J_\chi^\mathsf T\mathbf H_q\mathbf J_\chi
++ \sum_k (\mathbf g_q)_k\,\nabla^2_{\mathbf u}\chi_k.
+$$
+
+The final term is the geometric stiffness of the map and is included automatically.
+It can make the element Hessian indefinite even for a convex material, so use a projection mode that projects element Hessians to PD as needed, such as the default `ProjectionToPD::ProjectedNewton` (see [Newton's method](newtons_method.md)).
+Maps that couple a node's position to other nodes' DoFs, such as skinning or reduced coordinates, are not supported and are rejected; use `add_potential()` for those.
+
+
+## Indexable Summation Data
+
+In the previous Stable Neo-Hookean example, all elements share the same integration rule.
+Some applications need different weights, coordinates, or material data for each integration point of each element.
+SymX supports these inputs through `make_quadrature_scalar()`, `make_quadrature_vector()`, and `make_quadrature_matrix()`.
+Unlike `make_summation_vector(rule)`, these functions bind element-specific point data rather than a shared rule.
+
+Declare a fixed, positive number of points per element with `make_summation(count)`, or use the equivalent `add_for_each(count, definition)` wrapper for a nested loop style.
+Adapting the previous scalar potential definition:
+
+```cpp
+// Input stored outside the potential definition; populate before evaluation.
+std::vector<std::array<double, 4>> all_element_ip_rules; // {weight, xi0, xi1, xi2}
+
+// Inside the definition, using x, X, E, and nu as before:
+Scalar P = mws.add_for_each(n_integration_points_per_element,
+    [&]()
+    {
+        Vector ip = mws.make_quadrature_vector(all_element_ip_rules);
+        Scalar w = ip[0];
+        Vector xi({ip[1], ip[2], ip[3]});
+        Matrix Jx = fem_jacobian(data.element_type, x, xi);
+        Matrix JX = fem_jacobian(data.element_type, X, xi);
+        Matrix F = Jx * JX.inv();
+        return stable_neohookean_energy_density(F, E, nu) * JX.det() * w;
+    }
+);
+return P;
+```
+
+The array stores points in `[element][point][component]` order: point `p` of element `e` is entry `e * n_integration_points_per_element + p`.
+It must contain exactly `n_elements * n_integration_points_per_element` entries, with four components per entry in this example.
+All elements use the same point count, but their data may differ.
+Keep the input container alive for all evaluations.
+Indexable summation data is compatible with factorized FEM definitions.
+
+
+## Variable-Length Summation
+
+SymX also supports different summation lengths per element, for example to use integration rules with different numbers of points.
+Store the points consecutively as before, and pass an offsets array to `make_summation(offsets)` or `add_for_each(offsets, definition)` instead of a fixed count:
+
+```cpp
+// Input stored outside the potential definition, alongside all_element_ip_rules.
+std::vector<int64_t> offsets = {0, 2, 2, 5};
+// Populate all_element_ip_rules with five {weight, xi0, xi1, xi2} entries.
+
+// Inside the definition, using x, X, E, and nu as before:
+Scalar P = mws.add_for_each(offsets,
+    [&]()
+    {
+        Vector ip = mws.make_quadrature_vector(all_element_ip_rules);
+        Scalar w = ip[0];
+        Vector xi({ip[1], ip[2], ip[3]});
+        Matrix Jx = fem_jacobian(data.element_type, x, xi);
+        Matrix JX = fem_jacobian(data.element_type, X, xi);
+        Matrix F = Jx * JX.inv();
+        return stable_neohookean_energy_density(F, E, nu) * JX.det() * w;
+    }
+);
+return P;
+```
+
+Element `e` uses point entries in `[offsets[e], offsets[e + 1])`.
+Here the three elements have two, zero, and three points, respectively.
+An empty range contributes zero energy, gradient, and Hessian without evaluating the point kernel.
+Use `n_elements + 1` nondecreasing offsets starting at zero, with `offsets.back()` matching the number of point entries in each input array.
+
+Keep the offsets and point containers alive for all evaluations.
+You can update or resize them between evaluations without recompiling, provided the arrays remain consistent; do not modify them during evaluation.
+Variable-length summations also work with factorized FEM definitions.
+They support scalar `float` and `double` compiled kernels; SIMD compiled types are not supported.

docs/source/examples.md

@@ -84,6 +84,7 @@ Useful for understanding performance and the trade-offs between assembler back-e
 
 ---
 
+(rubber-extrusion-example)=
 ## Rubber Extrusion
 
 ![Rubber extrusion](_static/rubber_extrusion.jpg)
@@ -101,6 +102,10 @@ Key things to observe:
 
 Showcases how to use the high-level `NewtonsMethod` to solve FEM problems compactly and with minimal involvement from the user.
 
+The `RubberExtrusionMode` argument selects `StandardFixedSummation` for a shared rule, `StandardIndexedSummation` for indexed integration-point data, or `FactorizedFEM` for a split energy definition.
+The call in [`examples_main.cpp`](../../examples/examples_main.cpp) selects `FactorizedFEM`; change this argument to compare the definitions.
+See [Advanced FEM and Integration](advanced_integration.md) for the corresponding APIs and for variable-length summations.
+
 ---
 
 ## Inverse Cantilever

docs/source/fem_integration.md

@@ -99,6 +99,7 @@ return
 As described in [Symbol-Data Maps](symbol_data_maps.md), `MappedWorkspace` provides _summation_ facilities that allows for easy and efficient implementation of numerical integration.
 Note that SymX's summation is also _symbolic_ and **not** unrolled, which would explode the number of operations.
 Differentiation is efficiently handled as a summation.
+For element-specific integration-point data, variable point counts, and factorized FEM energies with mapped node positions, see [Advanced FEM and Integration](advanced_integration.md).
 
 ## FEM Integrator
 Finally, the whole summation logic is wrapped in `fem_integrator` for convenience:
@@ -118,6 +119,7 @@ Scalar fem_integrator(MappedWorkspace<double>& mws, const FEM_Element& element,
 }
 ```
 With this, the user only needs to indicate the element type and what to do with the integration point.
+The [rubber extrusion showcase](rubber-extrusion-example) includes a complete implementation using this integrator alongside indexed-quadrature and factorized FEM alternatives.
 
 
 ## Example

docs/source/global_optimization.md

@@ -120,3 +120,21 @@ G->apply_dof_increment(du.data()); // in-place update: x ← x + Δx
 ```
 
 These operations are used internally at each Newton iteration, and they are accessible to the user.
+
+## Energy and gradient breakdown
+
+Both `NewtonsMethod` and `FirstOrderOptimizer` expose their compiled evaluator through
+`get_compiled()`. To inspect individual energy contributions at the current DoFs:
+
+```cpp
+auto compiled = optimizer->get_compiled();
+auto energies = compiled->evaluate_P_components();
+auto components = compiled->evaluate_P__dP_du_components();
+
+double energy = components.energies.at("new_potential");
+const Eigen::VectorXd& gradient = components.gradients.at("new_potential");
+```
+
+Keys are the unique names passed to `add_potential()` or `add_potential_fem()`.
+Summing the component energies and gradients gives the results of
+`evaluate_P()` and `evaluate_P__dP_du()`.

docs/source/index.md

@@ -61,6 +61,7 @@ global_optimization
 fem_integration
 newtons_method
 first_order_optimization
+advanced_integration
 ```
 
 ```{toctree}

docs/source/newtons_method.md

@@ -21,11 +21,13 @@ Finally, the execution can occur.
 Because SymX has access to the DoF arrays (pointer and size), it can autonomously update their values in-place until the solution is reached.
 Multiple sets of degrees of freedom are allowed.
 
+For factorized FEM energies, mapped node positions, and integration data or point counts that vary by element, see [Advanced FEM and Integration](advanced_integration.md).
+
 ## Example
 We will build a 3D static elasticity solver using a FEM discretization for the Stable NeoHookean constitutive model and boundary conditions with the penalty method.
 The solver is compatible with the natively supported FEM element types SymX supports: `Tet3, Hex8, Tet10, Hex27`.
 Note that we use _Stable_ NeoHookean for the upcoming simulation examples for improved stability.
-A demo of solver can be found in [`rubber_extrusion.cpp`](../../examples/rubber_extrusion.cpp).
+The [rubber extrusion showcase](rubber-extrusion-example) provides a complete solver in [`rubber_extrusion.cpp`](../../examples/rubber_extrusion.cpp), with shared-rule, indexed-quadrature, and factorized FEM definitions.
 
 ```cpp
 struct Data