Tutorial: 3D Static Linear Elasticity with Affine Reduced Order Model

Problem Overview

This problem formulates and solves a three-dimensional static linear elasticity boundary value problem on a heterogeneous rectangular block. The block is composed of two distinct material subregions with the same constitutive model (isotropic linear elasticity) but different stiffness magnitudes. The problem is parametric in the elastic constants — Young’s modulus \(E\) and Poisson’s ratio \(\nu\) — and admits an affine decomposition of the stiffness operator, making it well-suited for Reduced Order Modeling (ROM) via Proper Orthogonal Decomposition (POD).

The mathematical treatment proceeds in the following stages:

Define the physical domain, boundary conditions, and material properties.
State the governing equations in strong and weak form.
Explain how the domain is partitioned into subregions and how their FE bases are constructed.
Derive the affine decomposition of the bilinear form, including the unit-coefficient assembly mechanism.
Describe the finite element discretization and natural boundary condition enforcement.
Formulate the Full Order Model (FOM).
Formulate the POD-based Reduced Order Model (ROM) with mean-shift correction.
Note the planned hyper-reduction extensions.

Link

Video tutorial

Project Architecture

The simulation workflow is organized through the following directory structure and files:

problem_1/
├─ domain.py           # Mesh generation, basis setup, boundary DOFs, and region partitioning
├─ bilinear_forms.py   # Unit-coefficient affine weak forms for λ and μ
├─ linear_forms.py     # Neumann traction right-hand-side term
├─ properties.py       # Regionwise Lamé parameter computation
├─ params.py           # Sobol sampling of (E, ν) for training and testing
├─ problem_def.py      # Main problem class registered with scikit-rom
└─ problem_1.ipynb     # Executes FOM generation, POD basis construction, ROM simulation, reporting, and VTK export

Physical Domain

Geometry

The computational domain is the rectangular parallelepiped:

\[\Omega = (0,\, L_x) \times (0,\, L_y) \times (0,\, L_z)\]

with physical dimensions \(L_x = 10\), \(L_y = 1\), \(L_z = 1\). The geometry resembles a slender beam-like block clamped at one end and loaded at the other.

Boundary Partitioning

The boundary \(\partial\Omega\) is partitioned into three disjoint parts:

Boundary	Location	Mechanical Condition
\(\Gamma_D\)	\(x = 0\) (left face)	Clamped: \(\mathbf{u} = \mathbf{0}\) (essential BC)
\(\Gamma_N\)	\(x = L_x = 10\) (right face)	Neumann: prescribed traction \(\mathbf{t}\) (natural BC)
\(\Gamma_\text{free}\)	top, bottom, front, back faces	Traction-free: \(\boldsymbol{\sigma}\mathbf{n} = \mathbf{0}\) (natural BC)

The essential (Dirichlet) condition on \(\Gamma_D\) is enforced explicitly by eliminating the corresponding degrees of freedom from the system. The natural conditions on \(\Gamma_N\) and \(\Gamma_\text{free}\) arise automatically from the variational formulation and require no additional enforcement beyond assembling the appropriate boundary integrals.

Material Subregions

The domain is partitioned into two non-overlapping subregions at the midpoint \(x = L_x/2 = 5\):

\[\Omega_1 = \{(x, y, z) \in \Omega : x \leq 5\}, \qquad \Omega_2 = \{(x, y, z) \in \Omega : x > 5\}\]

so that \(\Omega = \Omega_1 \cup \Omega_2\) and \(\Omega_1 \cap \Omega_2 = \emptyset\) (up to the shared interface at \(x=5\)). Region \(\Omega_1\) (the left half) is assigned a stiffness 50 times larger than Region \(\Omega_2\) (the right half), for all parameter values \((E, \nu)\).

Governing Equations

Kinematics: Linearized Strain Tensor

For small deformations, the linearized (symmetric) strain tensor is the symmetric part of the displacement gradient:

\[\boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\right)\]

Component-wise, for displacement \(\mathbf{u} = (u_1, u_2, u_3)^T\):

\[\varepsilon_{ij}(\mathbf{u}) = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right), \quad i,j \in \{1,2,3\}\]

The volumetric strain (trace of \(\boldsymbol{\varepsilon}\)) is:

\[\mathrm{tr}(\boldsymbol{\varepsilon}) = \varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33} = \nabla \cdot \mathbf{u}\]

Constitutive Law: Isotropic Linear Elasticity

The material obeys Hooke’s law for isotropic linear elastic solids. The Cauchy stress tensor \(\boldsymbol{\sigma}\) is related to the strain tensor by:

\[\boldsymbol{\sigma}(\mathbf{u}) = \lambda\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathbf{I} + 2\mu\,\boldsymbol{\varepsilon}(\mathbf{u})\]

where \(\mathbf{I}\) is the second-order identity tensor, and \(\lambda\), \(\mu\) are the Lamé parameters. The component-wise form is:

\[\sigma_{ij} = \lambda\,\varepsilon_{kk}\,\delta_{ij} + 2\mu\,\varepsilon_{ij}\]

where \(\delta_{ij}\) is the Kronecker delta and the repeated index \(k\) denotes summation (Einstein convention).

Lamé Parameters

Given Young’s modulus \(E > 0\) and Poisson’s ratio \(\nu \in (0, 0.5)\), the Lamé parameters are:

\[\lambda(E, \nu) = \frac{E\nu}{(1+\nu)(1-2\nu)}\]

\[\mu(E, \nu) = \frac{E}{2(1+\nu)}\]

Here \(\mu\) is the shear modulus and \(\lambda\) is the first Lamé parameter. Both are strictly positive for the admissible parameter range \(E \in [2, 10]\), \(\nu \in [0.2, 0.4]\).

Strong Form of the Boundary Value Problem

Find the displacement field \(\mathbf{u}: \Omega \to \mathbb{R}^3\) satisfying the equilibrium equation with zero body forces:

\[-\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) = \mathbf{0} \quad \text{in } \Omega\]

subject to the boundary conditions:

\[\mathbf{u} = \mathbf{0} \quad \text{on } \Gamma_D\]

\[\boldsymbol{\sigma}(\mathbf{u})\,\mathbf{n} = \mathbf{t} \quad \text{on } \Gamma_N\]

\[\boldsymbol{\sigma}(\mathbf{u})\,\mathbf{n} = \mathbf{0} \quad \text{on } \Gamma_\text{free}\]

where \(\mathbf{n}\) is the outward unit normal to \(\partial\Omega\).

Prescribed Traction

A uniform traction of magnitude \(0.01\) is applied in the negative \(y\)-direction on the right face \(\Gamma_N\):

\[\mathbf{t} = \begin{pmatrix} 0 \\ -0.01 \\ 0 \end{pmatrix}\]

This corresponds to a downward shear load applied at the free end of the block.

Weak Formulation

Function Spaces

Define the energy space for admissible displacements:

\[V = \left\{\mathbf{v} \in \left[H^1(\Omega)\right]^3 : \mathbf{v}\big|_{\Gamma_D} = \mathbf{0}\right\}\]

where \([H^1(\Omega)]^3\) denotes the Sobolev space of square-integrable vector fields with square-integrable first derivatives. This space encodes both the regularity required for the strain energy to be finite and the homogeneous Dirichlet condition.

Variational Problem

Multiply the equilibrium equation \(-\nabla \cdot \boldsymbol{\sigma} = \mathbf{0}\) by a test function \(\mathbf{v} \in V\), integrate over \(\Omega\), and apply the divergence theorem together with the boundary conditions. This yields the weak (variational) formulation:

Find \(\mathbf{u} \in V\) such that

\[a(\mathbf{u},\, \mathbf{v}) = \ell(\mathbf{v}) \qquad \forall\, \mathbf{v} \in V\]

Bilinear Form

The bilinear form \(a: V \times V \to \mathbb{R}\) represents the virtual work of internal stresses:

\[a(\mathbf{u}, \mathbf{v}) = \int_\Omega \boldsymbol{\sigma}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v})\, d\Omega\]

where the double contraction (Frobenius inner product) of two second-order tensors \(\mathbf{A}\) and \(\mathbf{B}\) is:

\[\mathbf{A} : \mathbf{B} = \sum_{i=1}^3 \sum_{j=1}^3 A_{ij}\, B_{ij}\]

Substituting the constitutive law and using the identity

\[\left(\lambda\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathbf{I}\right) : \boldsymbol{\varepsilon}(\mathbf{v}) = \lambda\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{v}))\]

the bilinear form expands to:

\[a(\mathbf{u}, \mathbf{v}) = \int_\Omega \left[\lambda\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{v})) + 2\mu\,\boldsymbol{\varepsilon}(\mathbf{u}):\boldsymbol{\varepsilon}(\mathbf{v})\right] d\Omega\]

This form is symmetric, continuous, and coercive on \(V\), guaranteeing a unique solution by the Lax–Milgram theorem.

Linear Form

The linear form \(\ell: V \to \mathbb{R}\) represents the virtual work of the external traction. Since \(\Gamma_\text{free}\) carries zero traction and \(\Gamma_D\) contributes nothing (test functions vanish there), only the right-face integral survives:

\[\ell(\mathbf{v}) = \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v}\, d\Gamma = -0.01 \int_{\Gamma_N} v_2\, d\Gamma\]

where \(v_2\) is the \(y\)-component of the test function.

Natural Boundary Condition Enforcement for Traction-Free Faces

The four remaining faces — top, bottom, front, and back (\(\Gamma_\text{free}\)) — carry zero traction. In the finite element method, these faces require no explicit enforcement: the condition \(\boldsymbol{\sigma}\,\mathbf{n} = \mathbf{0}\) is the natural boundary condition of the variational formulation, satisfied automatically when the corresponding boundary integral is omitted. Specifically:

A separate facet basis (fbasis_faces) is constructed over \(\{\)top, bottom, front, back\(\}\), and the corresponding DOF indices are identified.
No constraint is applied to these DOFs. The system is assembled without any traction contribution from these faces, which is equivalent to setting \(\mathbf{t} = \mathbf{0}\) there.
The traction-free condition thus drops out of the weak form at no additional cost.

Domain Partitioning and Region-Local FE Bases

Element-to-Location Mapping

To assemble the affine stiffness blocks separately over \(\Omega_1\) and \(\Omega_2\), the code must identify which elements belong to each subregion. This is done via the function element2location, which computes the physical \(x\)-coordinates of all quadrature/element nodes:

\[\text{loc}[e, g, 0] = x\text{-coordinate of quadrature point } g \text{ in element } e\]

An element belongs to Region 1 if all its nodal \(x\)-coordinates satisfy \(x \leq 5\) (the midpoint break):

\[\text{mask}_1[e] = \bigwedge_{g} \left(\text{loc}[e, g, 0] \leq \tfrac{L_x}{2}\right), \qquad \text{mask}_2[e] = \neg\,\text{mask}_1[e]\]

That is, an element is assigned to \(\Omega_1\) only if all of its nodes lie in the left half; otherwise it belongs to \(\Omega_2\).

Region-Local Basis Objects

Given the boolean element masks \(\text{mask}_i\), the function compute_basis_regions constructs a separate FE Basis object restricted to each subregion:

\[\texttt{basis\_regions} = \left\{\,\texttt{region\_1}: \mathcal{B}^{(1)},\quad \texttt{region\_2}: \mathcal{B}^{(2)}\,\right\}\]

where \(\mathcal{B}^{(i)}\) is the FE basis assembled over only the elements in \(\Omega_i\). Assembly over \(\mathcal{B}^{(i)}\) integrates only over the elements of the corresponding subregion, producing the parameter-independent stiffness blocks \(\mathbf{K}_q\) needed for the affine decomposition.

Material Heterogeneity

For a given parameter pair \(\boldsymbol{\mu} = (E, \nu)\), the region-wise Lamé parameters are:

Region 1 (\(\Omega_1\), stiffer left half):

\[\lambda_1(E,\nu) = 50\,\lambda(E,\nu), \qquad \mu_1(E,\nu) = 50\,\mu(E,\nu)\]

Region 2 (\(\Omega_2\), softer right half):

\[\lambda_2(E,\nu) = \lambda(E,\nu), \qquad \mu_2(E,\nu) = \mu(E,\nu)\]

The stiffness ratio between the two subregions is exactly \(50:1\) for all \((E, \nu)\) in the parameter space. For any region not explicitly assigned (fallback), the base parameters \((\lambda, \mu)\) are returned unchanged.

Affine Decomposition of the Bilinear Form

The affine decomposition is the key structural property that enables efficient Reduced Order Modeling. It separates parameter-dependent scalar coefficients from parameter-independent integral operators.

Unit-Coefficient Assembly Mechanism

The four stiffness blocks \(\mathbf{K}_q\) are assembled without any parameter values: when the assembly routine is called without specifying \(E\) or \(\nu\), the bilinear forms fall back to unit Lamé coefficients \(\lambda = 1\), \(\mu = 1\). This yields four purely geometric, parameter-free matrices:

\[(\mathbf{K}_\lambda^{(i)})_{jk} = \int_{\Omega_i} \mathrm{tr}(\boldsymbol{\varepsilon}(\boldsymbol{\phi}_k))\,\mathrm{tr}(\boldsymbol{\varepsilon}(\boldsymbol{\phi}_j))\, d\Omega, \quad i = 1, 2\]

\[(\mathbf{K}_\mu^{(i)})_{jk} = \int_{\Omega_i} 2\,\boldsymbol{\varepsilon}(\boldsymbol{\phi}_k):\boldsymbol{\varepsilon}(\boldsymbol{\phi}_j)\, d\Omega, \quad i = 1, 2\]

where \(\{\boldsymbol{\phi}_j\}_{j=1}^N\) are the vector FE basis functions. These four matrices are assembled once during the offline stage (at iteration counter cur_itr == 0) and cached. All subsequent queries reuse these cached blocks without re-assembly.

The actual material stiffness for a given \((E, \nu)\) is then recovered by scaling these unit blocks with the appropriate Lamé coefficients:

\[\mathbf{K}(E,\nu) = \theta_1\,\mathbf{K}_\lambda^{(1)} + \theta_2\,\mathbf{K}_\mu^{(1)} + \theta_3\,\mathbf{K}_\lambda^{(2)} + \theta_4\,\mathbf{K}_\mu^{(2)}\]

This unit-coefficient assembly is precisely the mechanism that makes the affine decomposition exact and exact at no approximation cost.

Splitting Over Subregions

Since \(\Omega = \Omega_1 \cup \Omega_2\), the bilinear form splits as:

\[a(\mathbf{u}, \mathbf{v};\, E, \nu) = a^{(1)}(\mathbf{u}, \mathbf{v};\, E, \nu) + a^{(2)}(\mathbf{u}, \mathbf{v};\, E, \nu)\]

where the contribution from subregion \(i\) is:

\[a^{(i)}(\mathbf{u}, \mathbf{v};\, E, \nu) = \int_{\Omega_i} \left[\lambda_i\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{v})) + 2\mu_i\,\boldsymbol{\varepsilon}(\mathbf{u}):\boldsymbol{\varepsilon}(\mathbf{v})\right] d\Omega\]

Parameter-Independent Sub-Forms

Define four parameter-independent bilinear sub-forms (assembled with unit coefficients as described above):

\[\mathcal{A}_\lambda^{(i)}(\mathbf{u}, \mathbf{v}) = \int_{\Omega_i} \mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{u}))\,\mathrm{tr}(\boldsymbol{\varepsilon}(\mathbf{v}))\, d\Omega, \quad i = 1, 2\]

\[\mathcal{A}_\mu^{(i)}(\mathbf{u}, \mathbf{v}) = \int_{\Omega_i} 2\,\boldsymbol{\varepsilon}(\mathbf{u}):\boldsymbol{\varepsilon}(\mathbf{v})\, d\Omega, \quad i = 1, 2\]

These integrals depend only on the geometry and mesh, not on \(E\) or \(\nu\).

Affine Expansion

Substituting the region-wise Lamé parameters, the full bilinear form admits the following \(Q = 4\) term affine expansion:

\[a(\mathbf{u}, \mathbf{v};\, E, \nu) = \underbrace{50\lambda(E,\nu)}_{\theta_1}\,\mathcal{A}_\lambda^{(1)} + \underbrace{50\mu(E,\nu)}_{\theta_2}\,\mathcal{A}_\mu^{(1)} + \underbrace{\lambda(E,\nu)}_{\theta_3}\,\mathcal{A}_\lambda^{(2)} + \underbrace{\mu(E,\nu)}_{\theta_4}\,\mathcal{A}_\mu^{(2)}\]

In compact notation, using \(\boldsymbol{\mu} = (E, \nu)\):

\[a(\mathbf{u}, \mathbf{v};\, \boldsymbol{\mu}) = \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathcal{A}_q(\mathbf{u}, \mathbf{v}), \qquad Q = 4\]

The four parameter-dependent coefficient functions and their corresponding sub-forms are:

Index \(q\)	Coefficient \(\theta_q(E, \nu)\)	Sub-form \(\mathcal{A}_q\)	Subregion
1	\(50\,\lambda(E,\nu)\)	\(\mathcal{A}_\lambda^{(1)}\) (\(\lambda\)-term)	\(\Omega_1\)
2	\(50\,\mu(E,\nu)\)	\(\mathcal{A}_\mu^{(1)}\) (\(\mu\)-term)	\(\Omega_1\)
3	\(\lambda(E,\nu)\)	\(\mathcal{A}_\lambda^{(2)}\) (\(\lambda\)-term)	\(\Omega_2\)
4	\(\mu(E,\nu)\)	\(\mathcal{A}_\mu^{(2)}\) (\(\mu\)-term)	\(\Omega_2\)

The linear form \(\ell(\mathbf{v})\) does not depend on the parameters \((E, \nu)\) and requires no affine decomposition.

Finite Element Discretization

Element Type

The domain \(\Omega\) is discretized with a structured mesh of trilinear hexahedral elements (Q1 elements). Each element has 8 nodes located at its corners. The isoparametric shape functions on the reference cube \(\hat{e} = [-1,1]^3\) are:

\[\hat{N}_I(\xi, \eta, \zeta) = \frac{1}{8}\left(1 + \xi_I\,\xi\right)\left(1 + \eta_I\,\eta\right)\left(1 + \zeta_I\,\zeta\right), \quad I = 1, \ldots, 8\]

where \((\xi_I, \eta_I, \zeta_I) \in \{-1,+1\}^3\) are the reference coordinates of the \(I\)-th corner node.

The displacement field within each element is approximated using a vector-valued basis — each nodal shape function multiplies each of the three displacement components independently:

\[\mathbf{u}^h(\mathbf{x}) = \sum_{I=1}^{8} \hat{N}_I\!\left(\boldsymbol{\xi}(\mathbf{x})\right)\,\mathbf{u}_I, \qquad \mathbf{u}_I \in \mathbb{R}^3\]

where \(\boldsymbol{\xi}(\mathbf{x})\) denotes the inverse isoparametric mapping from physical to reference coordinates.

The symmetric gradient entering the bilinear forms is computed on each element as:

\[\boldsymbol{\varepsilon}(\mathbf{u}^h)\big|_e = \frac{1}{2}\left(\nabla \mathbf{u}^h + (\nabla \mathbf{u}^h)^T\right)\big|_e\]

using the Jacobian of the isoparametric mapping to transform reference-element derivatives to physical derivatives.

Mesh Parameters

The mesh is generated with 8 divisions per unit length (refinement factor factor = 8), giving:

\[n_x = 8 \times L_x = 80, \quad n_y = 8 \times L_y = 8, \quad n_z = 8 \times L_z = 8\]

Quantity	Value
Block dimensions \((L_x \times L_y \times L_z)\)	\(10 \times 1 \times 1\)
Elements in \(x\)-direction	80
Elements in \(y\)-direction	8
Elements in \(z\)-direction	8
Total hexahedral elements	\(80 \times 8 \times 8 = 5{,}120\)
Nodes in \(x\)-direction	81
Nodes in \(y\)-direction	9
Nodes in \(z\)-direction	9
Total nodes	\(81 \times 9 \times 9 = 6{,}561\)
DOFs per node	3
Total degrees of freedom \(N\)	\(6{,}561 \times 3 = 19{,}683\)

Discrete System

After FE assembly, the variational problem reduces to the parametric linear algebraic system:

\[\mathbf{K}(\boldsymbol{\mu})\,\mathbf{u} = \mathbf{f}\]

where:

\(\mathbf{K}(\boldsymbol{\mu}) \in \mathbb{R}^{N \times N}\) is the global sparse stiffness matrix,
\(\mathbf{f} \in \mathbb{R}^N\) is the global load vector (assembled from the Neumann traction over \(\Gamma_N\)),
\(\mathbf{u} \in \mathbb{R}^N\) is the global displacement vector.

The affine structure is preserved exactly at the discrete level:

\[\mathbf{K}(\boldsymbol{\mu}) = \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathbf{K}_q, \qquad Q = 4\]

where each \(\mathbf{K}_q \in \mathbb{R}^{N \times N}\) is a sparse, symmetric, positive semi-definite matrix assembled over the corresponding subregion with unit Lamé coefficients, as described in the affine decomposition section.

The load vector is:

\[f_j = -0.01 \int_{\Gamma_N} \phi_j^{(2)}\, d\Gamma\]

where \(\phi_j^{(2)}\) denotes the \(y\)-component of the \(j\)-th vector shape function, integrated over the right-face facets.

Dirichlet Boundary Conditions

The homogeneous Dirichlet condition \(\mathbf{u} = \mathbf{0}\) on \(\Gamma_D\) is enforced by static condensation (elimination). The DOF indices on \(\Gamma_D\) are identified and the corresponding rows and columns are eliminated from the system. The condensed system involves only the free DOFs:

\[\widetilde{\mathbf{K}}(\boldsymbol{\mu})\,\tilde{\mathbf{u}} = \tilde{\mathbf{f}}\]

and is solved by a direct sparse solver. Since the Dirichlet value is zero (\(\mathbf{u}|_{\Gamma_D} = \mathbf{0}\)), no lift correction is needed on the right-hand side.

Parameter Space

The problem is parametrized by two material constants:

\[\boldsymbol{\mu} = (E,\, \nu) \in \mathcal{P} = [2,\, 10] \times [0.2,\, 0.4]\]

Sobol Sampling

Parameter samples are drawn using a Sobol quasi-random sequence, which is a low-discrepancy sequence that provides superior space-filling coverage of \(\mathcal{P}\) compared to random (Monte Carlo) sampling. The combined parameter array is arranged as training samples followed by testing samples, with boolean masks distinguishing the two sets:

Dataset	Number of Samples	Method
Training set \(\mathcal{P}_\text{train}\)	32	Sobol sequence
Testing set \(\mathcal{P}_\text{test}\)	32	Sobol sequence
Total	64	—

The training set is used in the offline phase to construct the reduced basis. The testing set is used in the online phase to evaluate ROM accuracy at unseen parameters.

Full Order Model (FOM)

The FOM solve for a given parameter \(\boldsymbol{\mu} = (E, \nu)\) proceeds as follows.

Step 1 — Assemble and Cache Affine Stiffness Blocks (Offline, Once)

At the first iteration (cur_itr == 0), load the domain data (mesh, basis, DOF indices, region bases), then assemble the four parameter-independent stiffness blocks over their respective region-local bases using unit Lamé coefficients. These blocks are stored in a nested dictionary K_list[region][material] and reused for all subsequent parameter queries without re-assembly:

\[\mathbf{K}_\lambda^{(i)} = \mathrm{asm}\!\left(\mathcal{A}_\lambda^{(i)},\; \mathcal{B}^{(i)}\right), \qquad \mathbf{K}_\mu^{(i)} = \mathrm{asm}\!\left(\mathcal{A}_\mu^{(i)},\; \mathcal{B}^{(i)}\right), \qquad i = 1, 2\]

Step 2 — Assemble and Cache Load Vector (Offline, Once)

Assemble the Neumann load vector by integrating the traction over the right-face facet basis:

\[f_j = \int_{\Gamma_N} \mathbf{t} \cdot \boldsymbol{\phi}_j\, d\Gamma = -0.01 \int_{\Gamma_N} \phi_j^{(2)}\, d\Gamma\]

This vector is also cached after the first assembly.

Step 3 — Combine Stiffness Matrix (Online, Per Query)

For each new \((E, \nu)\), compute the Lamé parameters and scale the cached unit blocks:

\[\mathbf{K}(E,\nu) = 50\lambda\,\mathbf{K}_\lambda^{(1)} + 50\mu\,\mathbf{K}_\mu^{(1)} + \lambda\,\mathbf{K}_\lambda^{(2)} + \mu\,\mathbf{K}_\mu^{(2)}\]

Step 4 — Apply Boundary Conditions and Solve

Initialize \(\mathbf{u} = \mathbf{0}\), set \(u[\Gamma_D] = 0\), condense out the \(\Gamma_D\) DOFs, and solve the condensed system by a direct sparse solver. The full displacement vector \(\mathbf{u}^h \in \mathbb{R}^N\) is returned as the FOM solution.

Reduced Order Model (ROM)

Offline Phase: Snapshot Collection and POD Basis Construction

Step 1 — Snapshot Collection

Solve the FOM for every training parameter sample:

\[\mathcal{S} = \left\{\mathbf{u}(\boldsymbol{\mu}^{(i)})\right\}_{i=1}^{N_\text{train}}, \qquad N_\text{train} = 32\]

Each snapshot \(\mathbf{u}(\boldsymbol{\mu}^{(i)}) \in \mathbb{R}^N\) is the full displacement field at training parameter \(\boldsymbol{\mu}^{(i)} = (E^{(i)}, \nu^{(i)})\).

Step 2 — Compute Mean Snapshot (Reference / Mean-Shift Vector)

\[\mathbf{T}_k = \frac{1}{N_\text{train}} \sum_{i=1}^{N_\text{train}} \mathbf{u}(\boldsymbol{\mu}^{(i)}) \in \mathbb{R}^N\]

This vector serves as the mean-shift (training reference solution). It is stored and used in both the offline projection step and the online right-hand side correction.

Step 3 — Form Mean-Centered Snapshot Matrix

\[\mathbf{S} = \left[\mathbf{u}(\boldsymbol{\mu}^{(1)}) - \mathbf{T}_k \;\Big|\; \mathbf{u}(\boldsymbol{\mu}^{(2)}) - \mathbf{T}_k \;\Big|\; \cdots \;\Big|\; \mathbf{u}(\boldsymbol{\mu}^{(N_\text{train})}) - \mathbf{T}_k\right] \in \mathbb{R}^{N \times N_\text{train}}\]

Step 4 — Singular Value Decomposition (POD)

Compute the thin SVD of the snapshot matrix:

\[\mathbf{S} = \mathbf{V}_\text{full}\,\boldsymbol{\Sigma}\,\mathbf{W}^T\]

where \(\mathbf{V}_\text{full} \in \mathbb{R}^{N \times N_\text{train}}\) contains the left singular vectors (POD modes), \(\boldsymbol{\Sigma} = \mathrm{diag}(\sigma_1, \ldots, \sigma_{N_\text{train}})\) with \(\sigma_1 \geq \sigma_2 \geq \cdots \geq 0\).

Retain the leading \(n \ll N\) modes to form the reduced basis matrix:

\[\mathbf{V} = [\mathbf{v}_1\;|\;\mathbf{v}_2\;|\;\cdots\;|\;\mathbf{v}_n] \in \mathbb{R}^{N \times n}\]

The POD modes are orthonormal: \(\mathbf{V}^T \mathbf{V} = \mathbf{I}_n\). The fraction of energy (variance) retained by the first \(n\) modes is:

\[\eta_n = \frac{\displaystyle\sum_{i=1}^{n} \sigma_i^2}{\displaystyle\sum_{i=1}^{N_\text{train}} \sigma_i^2}\]

A standard criterion is to choose \(n\) such that \(\eta_n \geq 1 - \varepsilon^2\) for a prescribed tolerance \(\varepsilon\).

Step 5 — Convert Stiffness Blocks to Sparse CSR Format

Before performing the matrix projections, each cached stiffness block is converted to Compressed Sparse Row (CSR) format:

\[\mathbf{K}_q^\text{CSR} = \mathrm{csr}(\mathbf{K}_q), \qquad q = 1, \ldots, Q\]

CSR format enables efficient sparse matrix–dense matrix products of the form \(\mathbf{V}^T \mathbf{K}_q^\text{CSR} \mathbf{V}\), which would otherwise be prohibitively expensive with unstructured sparse formats.

Step 6 — Project FOM Operators onto Reduced Basis

Precompute and store the following reduced-dimension quantities:

Projected stiffness blocks (\(n \times n\) dense matrices):

\[\mathbf{K}_{r,q} = \mathbf{V}^T \mathbf{K}_q\, \mathbf{V} \in \mathbb{R}^{n \times n}, \qquad q = 1, \ldots, Q\]

Projected mean-shift correction vectors (\(n\)-vectors):

\[\mathbf{c}_q = \mathbf{V}^T \mathbf{K}_q\, \mathbf{T}_k \in \mathbb{R}^n, \qquad q = 1, \ldots, Q\]

Projected load vector (\(n\)-vector):

\[\mathbf{f}_r = \mathbf{V}^T \mathbf{f} \in \mathbb{R}^n\]

All offline projections have cost \(O(N \cdot n)\) per block and are computed once during the offline stage (also triggered at cur_itr == 0).

Online Phase: ROM Solve

For a new query parameter \(\boldsymbol{\mu} = (E, \nu) \in \mathcal{P}\):

Step 1 — Compute Parameter-Dependent Coefficients

Evaluate the four affine coefficients using the Lamé formulas:

\[\theta_1 = 50\lambda(E,\nu), \quad \theta_2 = 50\mu(E,\nu), \quad \theta_3 = \lambda(E,\nu), \quad \theta_4 = \mu(E,\nu)\]

Step 2 — Assemble Reduced Stiffness Matrix

\[\mathbf{K}_r(\boldsymbol{\mu}) = \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathbf{K}_{r,q} \in \mathbb{R}^{n \times n}\]

Step 3 — Assemble Reduced Right-Hand Side with Mean-Shift Correction

\[\mathbf{g}(\boldsymbol{\mu}) = \mathbf{f}_r - \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathbf{c}_q \in \mathbb{R}^n\]

The subtracted term \(\sum_q \theta_q \mathbf{c}_q = \mathbf{V}^T \mathbf{K}(\boldsymbol{\mu})\,\mathbf{T}_k\) is the mean-shift correction: it accounts for the contribution of the mean snapshot \(\mathbf{T}_k\) to the residual after substituting the ROM ansatz.

Step 4 — Solve the Reduced System

\[\mathbf{K}_r(\boldsymbol{\mu})\,\mathbf{u}_r = \mathbf{g}(\boldsymbol{\mu})\]

This is a dense \(n \times n\) linear system solved by LU factorization at cost \(O(n^3)\). The result \(\mathbf{u}_r \in \mathbb{R}^n\) contains the reduced coordinates of the solution in the POD basis.

Step 5 — Reconstruct Full-Order Displacement (Framework Level)

The rom_solver method returns the reduced coordinates \(\mathbf{u}_r \in \mathbb{R}^n\). The full-order displacement reconstruction is performed externally by the calling framework:

\[\mathbf{u}_\text{ROM}(\boldsymbol{\mu}) = \mathbf{V}\,\mathbf{u}_r + \mathbf{T}_k \in \mathbb{R}^N\]

This expands the reduced solution back into the full \(N\)-dimensional space and adds the mean snapshot to account for the mean-shift. The separation between solving (inside rom_solver) and reconstructing (outside) reflects the library’s design: the framework handles post-processing uniformly for all problem types.

Derivation of the Mean-Shift Galerkin Projection

The ROM ansatz introduces the approximation:

\[\mathbf{u}(\boldsymbol{\mu}) \approx \mathbf{V}\,\mathbf{u}_r(\boldsymbol{\mu}) + \mathbf{T}_k\]

Substituting into the FOM system \(\mathbf{K}(\boldsymbol{\mu})\,\mathbf{u} = \mathbf{f}\):

\[\mathbf{K}(\boldsymbol{\mu})\left(\mathbf{V}\,\mathbf{u}_r + \mathbf{T}_k\right) = \mathbf{f}\]

\[\mathbf{K}(\boldsymbol{\mu})\,\mathbf{V}\,\mathbf{u}_r = \mathbf{f} - \mathbf{K}(\boldsymbol{\mu})\,\mathbf{T}_k\]

Applying the Galerkin projection by left-multiplying with \(\mathbf{V}^T\):

\[\underbrace{\mathbf{V}^T \mathbf{K}(\boldsymbol{\mu})\,\mathbf{V}}_{\mathbf{K}_r(\boldsymbol{\mu})}\,\mathbf{u}_r = \underbrace{\mathbf{V}^T \mathbf{f}}_{\mathbf{f}_r} - \underbrace{\mathbf{V}^T \mathbf{K}(\boldsymbol{\mu})\,\mathbf{T}_k}_{\sum_q \theta_q\,\mathbf{c}_q}\]

This is exactly the system \(\mathbf{K}_r(\boldsymbol{\mu})\,\mathbf{u}_r = \mathbf{g}(\boldsymbol{\mu})\) stated in the online phase. The affine structure is preserved throughout: since \(\mathbf{K}(\boldsymbol{\mu}) = \sum_q \theta_q \mathbf{K}_q\), each projected quantity inherits the same affine expansion:

\[\mathbf{K}_r(\boldsymbol{\mu}) = \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathbf{K}_{r,q}, \qquad \mathbf{V}^T \mathbf{K}(\boldsymbol{\mu})\,\mathbf{T}_k = \sum_{q=1}^{Q} \theta_q(\boldsymbol{\mu})\,\mathbf{c}_q\]

All terms on the right-hand side are precomputed offline, so the online cost is \(O(Qn^2)\) for assembly and \(O(n^3)\) for the solve — fully independent of the large FOM dimension \(N\).

Planned Hyper-Reduction Extensions

Three hyper-reduction methods are defined as placeholder stubs in the problem class, with full implementations planned for future development:

Method	Full Name	Purpose
DEIM	Discrete Empirical Interpolation Method	Efficient approximation of non-affine or nonlinear terms by interpolation at selected indices
ECSW	Energy-Conserving Sampling and Weighting	Reduces the cost of internal force assembly while preserving the energy structure
ECM	Empirical Cubature Method	Constructs a sparse cubature rule to replace full numerical integration in the ROM

These methods target the hyper-reduction stage: further reducing the online computational cost from \(O(Nn)\) (reconstruction) to \(O(pn)\) where \(p \ll N\) is the number of selected sample points or elements. The current implementation covers only the standard POD-Galerkin ROM without hyper-reduction.

Computational Cost Summary

Phase	FOM	ROM Offline	ROM Online
Domain & mesh setup	Once	Shared	—
Block assembly (\(\mathbf{K}_q\), \(\mathbf{f}\))	Once, cached	Shared	—
FOM solve per parameter	\(O(N^{1.5})\) sparse direct	\(N_\text{train}\) times	—
SVD / POD basis	—	\(O(N \cdot N_\text{train}^2)\)	—
Operator projection	—	\(O(Q \cdot N \cdot n)\)	—
Reduced system assembly	—	—	\(O(Qn^2)\)
Reduced solve	—	—	\(O(n^3)\)
Reconstruction	—	—	\(O(Nn)\)

Since \(n \ll N\) and \(Q = 4\), each online ROM evaluation is several orders of magnitude faster than a FOM solve while providing an accurate approximation of the full displacement field.

Summary of the Full Workflow

Offline stage (performed once):

Build the structured hexahedral mesh: \(80 \times 8 \times 8\) elements on \(\Omega = [0,10] \times [0,1] \times [0,1]\).
Partition elements into \(\Omega_1\) (\(x \leq 5\)) and \(\Omega_2\) (\(x > 5\)) using element2location; construct region-local FE bases via compute_basis_regions.
Assemble four parameter-independent unit-coefficient stiffness blocks \(\mathbf{K}_q\) over the region bases, and the load vector \(\mathbf{f}\) over the right-face facet basis. Cache all blocks.
Solve the FOM for \(N_\text{train} = 32\) Sobol-sampled training parameters \((E^{(i)}, \nu^{(i)}) \in [2,10] \times [0.2,0.4]\).
Compute the mean snapshot \(\mathbf{T}_k\) and the mean-centered snapshot matrix \(\mathbf{S}\).
Compute the thin SVD of \(\mathbf{S}\); select \(n\) POD modes as columns of \(\mathbf{V}\) based on the energy criterion \(\eta_n \geq 1 - \varepsilon^2\).
Convert stiffness blocks to CSR format; precompute and store \(\{\mathbf{K}_{r,q}\}\), \(\{\mathbf{c}_q\}\), and \(\mathbf{f}_r\).

Online stage (fast evaluation at any new \(\boldsymbol{\mu}\)):

Compute \(\theta_q(E,\nu)\) for \(q = 1,\ldots,4\).
Assemble \(\mathbf{K}_r(\boldsymbol{\mu}) = \sum_q \theta_q \mathbf{K}_{r,q}\) and \(\mathbf{g}(\boldsymbol{\mu}) = \mathbf{f}_r - \sum_q \theta_q \mathbf{c}_q\).
Solve the \(n \times n\) system \(\mathbf{K}_r \mathbf{u}_r = \mathbf{g}\) by LU factorization.
Reconstruct the full displacement: \(\mathbf{u}_\text{ROM} = \mathbf{V}\mathbf{u}_r + \mathbf{T}_k \in \mathbb{R}^N\).