Solid dynamics model#

The examples include a lightweight two-dimensional solid-dynamics solver written in pure Python with NumPy and SciPy. The model is intended as a transparent reference problem for projection-based ROMs, hyper-reduction, surrogates, UQ, and inverse workflows. It is not intended to replace a production structural-mechanics code.

The core Neo-Hookean solver, the common material-model facade, and the romtools QoI wrapper live in examples/models. The user-facing examples live in examples/nonlinear_solid_dynamics.

Governing equations#

After finite-element discretization the displacement satisfies

\[M \ddot{u} + C \dot{u} + f_{\mathrm{int}}(u;\mu) = f_{\mathrm{ext}}(t;\mu).\]

Both displacement and velocity are retained as state variables. The mechanics operators are exposed independently of the time integrator through methods such as mass_matrix, internal_force, tangent_stiffness, and acceleration. Element-level internal-force evaluation is also available for future hyper-reduction examples.

Spatial discretization and materials#

The implementation uses a structured rectangular mesh of bilinear Q4 elements, two-by-two Gauss quadrature, and plane-strain kinematics. Two constitutive models are available through examples/models/solid_dynamics.py.

Neo-Hookean#

The nonlinear model uses a total-Lagrangian formulation and the compressible Neo-Hookean energy

\[W(C,J) = \frac{\lambda}{4}(J^2-1) - \left(\frac{\lambda}{2}+\mu\right)\log J + \frac{\mu}{2}(I_C-3),\]

where lambda and mu are the Lamé constants. The corresponding first Piola stress is

\[P = \mu F + \left[\frac{\lambda}{2}(J^2-1)-\mu\right]F^{-T}.\]

The consistent material tangent is assembled analytically.

Linear elasticity#

A small-strain isotropic linear-elastic option is also available. It uses

\[\varepsilon = \frac{1}{2}\left(\nabla u + \nabla u^T\right),\]
\[\sigma = \lambda\,\mathrm{tr}(\varepsilon)I + 2\mu\varepsilon.\]

The resulting internal force is linear in displacement and the tangent stiffness is constant. The linear material reuses the same mesh, mass matrix, boundary-condition handling, and explicit/implicit time integrators as the Neo-Hookean material. The beam and doubly clamped examples accept material_model="neo_hookean" or material_model="linear".

Both consistent and row-sum lumped mass matrices are supported.

Time integration#

solve_explicit uses velocity Verlet. solve_implicit uses Newmark with the average-acceleration defaults

\[\gamma=\frac{1}{2}, \qquad \beta=\frac{1}{4}.\]

The implicit formulation deliberately uses velocity at the new time step as the Newton unknown. For a trial v_{n+1},

\[a_{n+1} = \frac{v_{n+1}-v_n}{\gamma\Delta t} - \frac{1-\gamma}{\gamma}a_n,\]

and u_{n+1} follows from the Newmark displacement relation. The residual is

\[R(v_{n+1}) = M a_{n+1} + C v_{n+1} + f_{\mathrm{int}}(u_{n+1}) - f_{\mathrm{ext}}(t_{n+1}),\]

with Jacobian

\[\frac{\partial R}{\partial v} = \frac{1}{\gamma\Delta t}M + C + \frac{\beta\Delta t}{\gamma}K_{\mathrm{tan}}.\]

Newton iterations use a simple backtracking line search. For the linear material the same formulation is retained, although the residual is affine and one exact Newton correction is sufficient up to numerical roundoff.

Examples#

Cantilever beam#

examples/nonlinear_solid_dynamics/beam.py clamps the left end of a rectangular beam and applies a smooth transient transverse resultant on the right end. The example plots the vertical tip-displacement history and stores full displacement/velocity trajectories in memory for direct use as ROM snapshots. Its run function accepts either constitutive model.

Run the default Neo-Hookean case with

python examples/nonlinear_solid_dynamics/beam.py

Doubly clamped Gaussian perturbation#

examples/nonlinear_solid_dynamics/clamped_gaussian.py fixes both ends and starts from a smooth transverse Gaussian displacement centered in the domain,

\[u_y(X,0) = A\exp\left[-\frac{(X-X_c)^2}{2\sigma^2}\right].\]

The amplitude, width, and center are explicit parameters of the initialization. The same example can use either the Neo-Hookean or linear-elastic material.

The linear transverse problem is not the scalar two-way wave equation. A transverse disturbance of a clamped beam excites flexural/elastic modes; in a slender-beam approximation the corresponding Euler–Bernoulli equation contains a fourth spatial derivative and is dispersive. In the full 2-D linear solid, the governing equations are the vector elastodynamic equations and support longitudinal and shear waves.

Longitudinal two-way wave check#

examples/nonlinear_solid_dynamics/linear_wave.py provides a separate configuration that does reduce to the scalar two-way wave equation. It uses the linear material, constrains all transverse DOFs, fixes the axial displacement at both ends, and initializes an axial Gaussian perturbation that is uniform through the height. The continuum reduction is

\[u_{tt} = c_p^2 u_{xx},\]

with the plane-strain longitudinal wave speed

\[c_p = \sqrt{\frac{\lambda+2\mu}{\rho}}.\]

For zero initial velocity and the Gaussian initial displacement

\[g(x)=A\exp\left[-\frac{(x-x_c)^2}{2\sigma^2}\right],\]

D’Alembert’s solution on the line is

\[u(x,t)=\frac{1}{2}g(x-c_pt)+\frac{1}{2}g(x+c_pt).\]

The example compares the finite-element solution directly with this analytic solution before either pulse reaches a clamped boundary. The initial Gaussian is chosen sufficiently far from the ends that its boundary value is negligible over the comparison interval.

The script also performs a coupled spatial/temporal refinement study with dt proportional to dx/c_p. It reports the relative discrete \(L^2\) displacement error and observed convergence rate for successively refined Q4 meshes. The corresponding CI test requires the error to decrease under refinement and checks for the expected near-second-order trend without using a platform-sensitive exact rate tolerance.

Run the analytic comparison and convergence study with

python examples/nonlinear_solid_dynamics/linear_wave.py

Published nonlinear benchmark check#

examples/nonlinear_solid_dynamics/benchmark_cantilever.py reproduces the large-deformation cantilever benchmark of Stickle et al. (2022). The published problem is a 4 m by 1 m plane-strain Neo-Hookean cantilever with

\[E=1\ \mathrm{MPa},\qquad \nu=0.3,\qquad \rho_0=1050\ \mathrm{kg/m^3},\qquad g=10\ \mathrm{m/s^2}.\]

The left end is clamped and gravity is applied instantaneously at t=0. The reference calculation uses Newmark with beta=0.25, gamma=0.5 and dt=1e-3 s. The paper uses an 8 by 2 Q8 mesh for the plotted benchmark; the romtools reference implementation intentionally uses Q4 elements.

For the 8 by 2 Q4 mesh, the present implementation gives a first lower-right vertical-displacement minimum of approximately -3.17 m at t=1.00 s. The published response reaches approximately -3.3 m near t=1.0--1.1 s. The agreement is sufficient for the intended lightweight implementation check, particularly given the difference in element order.

A short CI regression uses dt=1e-2 s through the first displacement minimum. The full example retains the published dt=1e-3 s and 3 s time horizon.

Verification scope#

The initial solver keeps verification deliberately focused:

  • zero internal force is checked in the undeformed configuration;

  • the Neo-Hookean analytical tangent is compared with a finite-difference directional derivative of the internal force;

  • the linear material is checked for force linearity and a constant tangent;

  • the longitudinal linear-wave configuration is checked for left/right symmetry and convergence to the analytic D’Alembert solution under mesh/time refinement;

  • a small velocity-primary Newmark solve checks constrained DOFs and nonlinear convergence; and

  • the first minimum of the Stickle cantilever response is used as an external nonlinear regression check.

These checks establish confidence in the example implementation without turning romtools into a general solid-mechanics verification suite.

QoI wrapper#

NonlinearSolidBeamQoiModel in examples/models/nonlinear_solid_dynamics_model.py exposes a parameterized cantilever through the standard romtools model protocol. The sample parameters are young_modulus, load_amplitude, and pulse_duration. The QoI is the vertical tip-displacement history; full state snapshots are saved alongside the QoI for later ROM construction. The wrapper constructor also accepts material_model="neo_hookean" or material_model="linear".

Implementation#

Neo-Hookean core solver:

  1"""Lightweight 2-D nonlinear solid dynamics reference model.
  2
  3This module implements a structured Q4 total-Lagrangian finite-element model
  4with a compressible Neo-Hookean material.  It is intentionally compact and is
  5meant for romtools examples, not production structural analysis.
  6"""
  7
  8from __future__ import annotations
  9
 10from dataclasses import dataclass
 11from typing import Callable, Optional
 12
 13import numpy as np
 14import scipy.linalg
 15
 16Array = np.ndarray
 17
 18
 19@dataclass(frozen=True)
 20class Material:
 21    """Compressible Neo-Hookean material in plane strain."""
 22
 23    young_modulus: float
 24    poisson_ratio: float
 25    density: float
 26
 27    @property
 28    def lame_lambda(self) -> float:
 29        e = float(self.young_modulus)
 30        nu = float(self.poisson_ratio)
 31        return e * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
 32
 33    @property
 34    def lame_mu(self) -> float:
 35        e = float(self.young_modulus)
 36        nu = float(self.poisson_ratio)
 37        return e / (2.0 * (1.0 + nu))
 38
 39
 40@dataclass
 41class SolidState:
 42    displacement: Array
 43    velocity: Array
 44    acceleration: Array
 45    time: float = 0.0
 46
 47
 48class StructuredQ4Mesh:
 49    """Uniform rectangular Q4 mesh."""
 50
 51    def __init__(self, length: float, height: float, nx: int, ny: int):
 52        if nx < 1 or ny < 1:
 53            raise ValueError("nx and ny must both be positive")
 54        self.length = float(length)
 55        self.height = float(height)
 56        self.nx = int(nx)
 57        self.ny = int(ny)
 58
 59        xs = np.linspace(0.0, self.length, self.nx + 1)
 60        ys = np.linspace(0.0, self.height, self.ny + 1)
 61        self.coordinates = np.array([[x, y] for y in ys for x in xs], dtype=float)
 62
 63        elements = []
 64        stride = self.nx + 1
 65        for j in range(self.ny):
 66            for i in range(self.nx):
 67                n0 = j * stride + i
 68                n1 = n0 + 1
 69                n3 = n0 + stride
 70                n2 = n3 + 1
 71                elements.append([n0, n1, n2, n3])
 72        self.elements = np.asarray(elements, dtype=int)
 73
 74    @property
 75    def num_nodes(self) -> int:
 76        return self.coordinates.shape[0]
 77
 78    @property
 79    def num_elements(self) -> int:
 80        return self.elements.shape[0]
 81
 82    @property
 83    def ndof(self) -> int:
 84        return 2 * self.num_nodes
 85
 86    def nodes_on_x(self, x: float, atol: float = 1.0e-12) -> Array:
 87        return np.flatnonzero(np.isclose(self.coordinates[:, 0], x, atol=atol, rtol=0.0))
 88
 89    def node_nearest(self, x: float, y: float) -> int:
 90        point = np.array([x, y], dtype=float)
 91        return int(np.argmin(np.linalg.norm(self.coordinates - point[None, :], axis=1)))
 92
 93
 94def _shape_q4(xi: float, eta: float) -> tuple[Array, Array]:
 95    n = 0.25 * np.array(
 96        [
 97            (1.0 - xi) * (1.0 - eta),
 98            (1.0 + xi) * (1.0 - eta),
 99            (1.0 + xi) * (1.0 + eta),
100            (1.0 - xi) * (1.0 + eta),
101        ]
102    )
103    dndxi = 0.25 * np.array(
104        [
105            [-(1.0 - eta), -(1.0 - xi)],
106            [+(1.0 - eta), -(1.0 + xi)],
107            [+(1.0 + eta), +(1.0 + xi)],
108            [-(1.0 + eta), +(1.0 - xi)],
109        ]
110    )
111    return n, dndxi
112
113
114_GAUSS = 1.0 / np.sqrt(3.0)
115_GAUSS_POINTS = ((-_GAUSS, -_GAUSS), (_GAUSS, -_GAUSS), (_GAUSS, _GAUSS), (-_GAUSS, _GAUSS))
116
117
118class NonlinearSolid2D:
119    """Total-Lagrangian Q4 finite-element solid with plane-strain kinematics."""
120
121    def __init__(
122        self,
123        mesh: StructuredQ4Mesh,
124        material: Material,
125        constrained_dofs: Optional[Array] = None,
126        mass_type: str = "consistent",
127        rayleigh_alpha: float = 0.0,
128        rayleigh_beta: float = 0.0,
129    ):
130        self.mesh = mesh
131        self.material = material
132        self.constrained_dofs = np.unique(
133            np.asarray([] if constrained_dofs is None else constrained_dofs, dtype=int)
134        )
135        all_dofs = np.arange(mesh.ndof, dtype=int)
136        self.free_dofs = np.setdiff1d(all_dofs, self.constrained_dofs, assume_unique=True)
137        if mass_type not in ("consistent", "lumped"):
138            raise ValueError("mass_type must be 'consistent' or 'lumped'")
139        self.mass_type = mass_type
140        self._quadrature_cache = [self._compute_element_quadrature(e) for e in range(mesh.num_elements)]
141        self._mass = self._assemble_mass(lumped=(mass_type == "lumped"))
142        self.rayleigh_alpha = float(rayleigh_alpha)
143        self.rayleigh_beta = float(rayleigh_beta)
144        self._initial_tangent = None
145
146    @staticmethod
147    def clamp_dofs(mesh: StructuredQ4Mesh, left: bool = True, right: bool = False) -> Array:
148        nodes = []
149        if left:
150            nodes.extend(mesh.nodes_on_x(0.0).tolist())
151        if right:
152            nodes.extend(mesh.nodes_on_x(mesh.length).tolist())
153        nodes = np.unique(np.asarray(nodes, dtype=int))
154        dofs = np.empty(2 * nodes.size, dtype=int)
155        dofs[0::2] = 2 * nodes
156        dofs[1::2] = 2 * nodes + 1
157        return np.sort(dofs)
158
159    @property
160    def ndof(self) -> int:
161        return self.mesh.ndof
162
163    def mass_matrix(self) -> Array:
164        return self._mass.copy()
165
166    def damping_matrix(self) -> Array:
167        c = self.rayleigh_alpha * self._mass
168        if self.rayleigh_beta != 0.0:
169            if self._initial_tangent is None:
170                self._initial_tangent = self.tangent_stiffness(np.zeros(self.ndof))
171            c = c + self.rayleigh_beta * self._initial_tangent
172        return c
173
174    def _element_dofs(self, nodes: Array) -> Array:
175        dofs = np.empty(2 * len(nodes), dtype=int)
176        dofs[0::2] = 2 * nodes
177        dofs[1::2] = 2 * nodes + 1
178        return dofs
179
180    def _compute_element_quadrature(self, element_id: int):
181        nodes = self.mesh.elements[element_id]
182        x = self.mesh.coordinates[nodes]
183        data = []
184        for xi, eta in _GAUSS_POINTS:
185            n, dndxi = _shape_q4(xi, eta)
186            jac = x.T @ dndxi
187            det_jac = np.linalg.det(jac)
188            if det_jac <= 0.0:
189                raise ValueError("reference element has non-positive Jacobian")
190            dndx = dndxi @ np.linalg.inv(jac)
191            data.append((n, dndx, det_jac))
192        return tuple(data)
193
194    def _element_quadrature(self, element_id: int):
195        return self._quadrature_cache[element_id]
196
197    def _assemble_mass(self, lumped: bool) -> Array:
198        m = np.zeros((self.ndof, self.ndof), dtype=float)
199        rho = self.material.density
200        for e, nodes in enumerate(self.mesh.elements):
201            me_scalar = np.zeros((4, 4), dtype=float)
202            for n, _dndx, dv in self._element_quadrature(e):
203                me_scalar += rho * np.outer(n, n) * dv
204            me = np.kron(me_scalar, np.eye(2))
205            dofs = self._element_dofs(nodes)
206            m[np.ix_(dofs, dofs)] += me
207        if lumped:
208            m = np.diag(np.sum(m, axis=1))
209        return m
210
211    def _first_piola(self, f2: Array) -> tuple[Array, float, Array, float]:
212        f = np.eye(3)
213        f[:2, :2] = f2
214        j = np.linalg.det(f)
215        if j <= 0.0:
216            raise ValueError("non-positive deformation Jacobian")
217        mu = self.material.lame_mu
218        lam = self.material.lame_lambda
219        finvt = np.linalg.inv(f).T
220        a = 0.5 * lam * (j * j - 1.0) - mu
221        p = mu * f + a * finvt
222        return p[:2, :2], j, finvt, a
223
224    def _first_piola_and_tangent(self, f2: Array) -> tuple[Array, Array]:
225        """Return in-plane P and dP/dF for the paper's compressible NH model."""
226        p, j, finvt, a = self._first_piola(f2)
227        mu = self.material.lame_mu
228        lam = self.material.lame_lambda
229        eye = np.eye(3)
230        amat = (
231            mu * np.einsum("ik,JL->iJkL", eye, eye)
232            + lam * j * j * np.einsum("iJ,kL->iJkL", finvt, finvt)
233            - a * np.einsum("kJ,iL->iJkL", finvt, finvt)
234        )
235        return p, amat[:2, :2, :2, :2]
236
237    def _element_force_and_tangent(self, element_id: int, u_local: Array) -> tuple[Array, Array]:
238        u_local = np.asarray(u_local, dtype=float).reshape(4, 2)
239        fe = np.zeros((4, 2), dtype=float)
240        ke4 = np.zeros((4, 2, 4, 2), dtype=float)
241        for _n, grad_n, dv in self._element_quadrature(element_id):
242            grad_u = u_local.T @ grad_n
243            f2 = np.eye(2) + grad_u
244            p, amat = self._first_piola_and_tangent(f2)
245            fe += (grad_n @ p.T) * dv
246            ke4 += np.einsum("aJ,iJkL,bL->aibk", grad_n, amat, grad_n) * dv
247        return fe.reshape(-1), ke4.reshape(8, 8)
248
249    def element_internal_force(self, element_id: int, u_local: Array) -> Array:
250        u_local = np.asarray(u_local, dtype=float).reshape(4, 2)
251        fe = np.zeros((4, 2), dtype=float)
252        for _n, grad_n, dv in self._element_quadrature(element_id):
253            grad_u = u_local.T @ grad_n
254            p, _j, _finvt, _a = self._first_piola(np.eye(2) + grad_u)
255            fe += (grad_n @ p.T) * dv
256        return fe.reshape(-1)
257
258    def element_tangent_stiffness(self, element_id: int, u_local: Array) -> Array:
259        return self._element_force_and_tangent(element_id, u_local)[1]
260
261    def internal_force_and_tangent(self, displacement: Array) -> tuple[Array, Array]:
262        u = np.asarray(displacement, dtype=float)
263        fint = np.zeros(self.ndof, dtype=float)
264        k = np.zeros((self.ndof, self.ndof), dtype=float)
265        for e, nodes in enumerate(self.mesh.elements):
266            dofs = self._element_dofs(nodes)
267            fe, ke = self._element_force_and_tangent(e, u[dofs])
268            fint[dofs] += fe
269            k[np.ix_(dofs, dofs)] += ke
270        return fint, k
271
272    def internal_force(self, displacement: Array) -> Array:
273        u = np.asarray(displacement, dtype=float)
274        fint = np.zeros(self.ndof, dtype=float)
275        for e, nodes in enumerate(self.mesh.elements):
276            dofs = self._element_dofs(nodes)
277            fint[dofs] += self.element_internal_force(e, u[dofs])
278        return fint
279
280    def tangent_stiffness(self, displacement: Array) -> Array:
281        return self.internal_force_and_tangent(displacement)[1]
282
283    def body_force(self, acceleration: Array) -> Array:
284        acc = np.asarray(acceleration, dtype=float)
285        if acc.shape != (2,):
286            raise ValueError("body acceleration must have shape (2,)")
287        nodal_acc = np.tile(acc, self.mesh.num_nodes)
288        return self._mass @ nodal_acc
289
290    def nodal_force(self, node: int, force: Array) -> Array:
291        out = np.zeros(self.ndof, dtype=float)
292        out[2*node:2*node+2] = np.asarray(force, dtype=float)
293        return out
294
295    def boundary_force_x(self, x: float, total_force: Array) -> Array:
296        """Distribute a resultant force consistently along a vertical boundary."""
297        nodes = self.mesh.nodes_on_x(x)
298        nodes = nodes[np.argsort(self.mesh.coordinates[nodes, 1])]
299        weights = np.zeros(nodes.size)
300        y = self.mesh.coordinates[nodes, 1]
301        for i in range(nodes.size - 1):
302            segment = y[i + 1] - y[i]
303            weights[i] += 0.5 * segment
304            weights[i + 1] += 0.5 * segment
305        weights /= np.sum(weights)
306        out = np.zeros(self.ndof)
307        force = np.asarray(total_force, dtype=float)
308        for node, weight in zip(nodes, weights):
309            out[2*node:2*node+2] += weight * force
310        return out
311
312    def acceleration(
313        self,
314        displacement: Array,
315        velocity: Array,
316        time: float,
317        external_force: Callable[[float], Array],
318    ) -> Array:
319        c = self.damping_matrix()
320        rhs = external_force(time) - self.internal_force(displacement) - c @ velocity
321        a = np.zeros(self.ndof, dtype=float)
322        ff = np.ix_(self.free_dofs, self.free_dofs)
323        a[self.free_dofs] = scipy.linalg.solve(self._mass[ff], rhs[self.free_dofs], assume_a="sym")
324        return a
325
326    def initial_state(
327        self,
328        displacement: Optional[Array] = None,
329        velocity: Optional[Array] = None,
330        external_force: Optional[Callable[[float], Array]] = None,
331        time: float = 0.0,
332    ) -> SolidState:
333        u = np.zeros(self.ndof) if displacement is None else np.array(displacement, dtype=float, copy=True)
334        v = np.zeros(self.ndof) if velocity is None else np.array(velocity, dtype=float, copy=True)
335        u[self.constrained_dofs] = 0.0
336        v[self.constrained_dofs] = 0.0
337        if external_force is None:
338            external_force = lambda _t: np.zeros(self.ndof)
339        a = self.acceleration(u, v, time, external_force)
340        return SolidState(u, v, a, float(time))
341
342    def solve_implicit(
343        self,
344        state: SolidState,
345        dt: float,
346        num_steps: int,
347        external_force: Callable[[float], Array],
348        beta: float = 0.25,
349        gamma: float = 0.5,
350        newton_tolerance: float = 1.0e-9,
351        max_newton_iterations: int = 25,
352        line_search: bool = True,
353        snapshot_stride: int = 1,
354    ) -> tuple[Array, Array, Array]:
355        """Integrate with velocity-primary Newmark and Newton iteration."""
356        u = np.array(state.displacement, copy=True)
357        v = np.array(state.velocity, copy=True)
358        a = np.array(state.acceleration, copy=True)
359        t = float(state.time)
360        cmat = self.damping_matrix()
361        ff = np.ix_(self.free_dofs, self.free_dofs)
362        mff = self._mass[ff]
363        cff = cmat[ff]
364        cu = beta * dt / gamma
365
366        times = [t]
367        displacements = [u.copy()]
368        velocities = [v.copy()]
369
370        for istep in range(1, num_steps + 1):
371            tnext = t + dt
372            v_old = v.copy()
373            u_old = u.copy()
374            a_old = a.copy()
375            v_trial = v_old + dt * a_old
376            v_trial[self.constrained_dofs] = 0.0
377
378            def kinematics(v_candidate):
379                anew = (
380                    (v_candidate - v_old) / (gamma * dt)
381                    - ((1.0 - gamma) / gamma) * a_old
382                )
383                unew = (
384                    u_old
385                    + dt * v_old
386                    + dt * dt * ((0.5 - beta) * a_old + beta * anew)
387                )
388                unew[self.constrained_dofs] = 0.0
389                anew[self.constrained_dofs] = 0.0
390                return unew, anew
391
392            converged = False
393            fext = external_force(tnext)
394            for _iteration in range(max_newton_iterations):
395                u_trial, a_trial = kinematics(v_trial)
396                fint, kt = self.internal_force_and_tangent(u_trial)
397                inertial = self._mass @ a_trial
398                damping = cmat @ v_trial
399                residual = inertial + damping + fint - fext
400                rfree = residual[self.free_dofs]
401                scale = max(
402                    np.linalg.norm(fext[self.free_dofs]),
403                    np.linalg.norm(inertial[self.free_dofs]),
404                    np.linalg.norm(fint[self.free_dofs]),
405                    1.0,
406                )
407                rnorm = np.linalg.norm(rfree)
408                if rnorm <= newton_tolerance * scale:
409                    converged = True
410                    break
411
412                jac = mff / (gamma * dt) + cff + cu * kt[ff]
413                dv = scipy.linalg.solve(jac, -rfree, assume_a="sym")
414
415                alpha = 1.0
416                if line_search:
417                    base = rnorm
418                    for _ in range(10):
419                        candidate = v_trial.copy()
420                        candidate[self.free_dofs] += alpha * dv
421                        candidate[self.constrained_dofs] = 0.0
422                        uc, ac = kinematics(candidate)
423                        rc = (
424                            self._mass @ ac
425                            + cmat @ candidate
426                            + self.internal_force(uc)
427                            - fext
428                        )[self.free_dofs]
429                        if np.linalg.norm(rc) < base:
430                            break
431                        alpha *= 0.5
432                v_trial[self.free_dofs] += alpha * dv
433                v_trial[self.constrained_dofs] = 0.0
434
435            if not converged:
436                raise RuntimeError(f"Newton failed to converge at step {istep}, t={tnext:.6e}")
437
438            u, a = kinematics(v_trial)
439            v = v_trial
440            t = tnext
441            if istep % snapshot_stride == 0 or istep == num_steps:
442                times.append(t)
443                displacements.append(u.copy())
444                velocities.append(v.copy())
445
446        return np.asarray(times), np.asarray(displacements), np.asarray(velocities)
447
448    def solve_explicit(
449        self,
450        state: SolidState,
451        dt: float,
452        num_steps: int,
453        external_force: Callable[[float], Array],
454        snapshot_stride: int = 1,
455    ) -> tuple[Array, Array, Array]:
456        """Integrate with velocity Verlet."""
457        u = np.array(state.displacement, copy=True)
458        v = np.array(state.velocity, copy=True)
459        a = np.array(state.acceleration, copy=True)
460        t = float(state.time)
461        times = [t]
462        displacements = [u.copy()]
463        velocities = [v.copy()]
464
465        for istep in range(1, num_steps + 1):
466            vhalf = v + 0.5 * dt * a
467            unew = u + dt * vhalf
468            unew[self.constrained_dofs] = 0.0
469            anew = self.acceleration(unew, vhalf, t + dt, external_force)
470            vnew = vhalf + 0.5 * dt * anew
471            vnew[self.constrained_dofs] = 0.0
472            u, v, a = unew, vnew, anew
473            t += dt
474            if istep % snapshot_stride == 0 or istep == num_steps:
475                times.append(t)
476                displacements.append(u.copy())
477                velocities.append(v.copy())
478        return np.asarray(times), np.asarray(displacements), np.asarray(velocities)
479
480    def gaussian_transverse_displacement(
481        self, amplitude: float, width: float, center: Optional[float] = None
482    ) -> Array:
483        if center is None:
484            center = 0.5 * self.mesh.length
485        x = self.mesh.coordinates[:, 0]
486        profile = float(amplitude) * np.exp(-0.5 * ((x - center) / float(width)) ** 2)
487        u = np.zeros(self.ndof)
488        u[1::2] = profile
489        u[self.constrained_dofs] = 0.0
490        return u
491
492    def gaussian_transverse_velocity(
493        self, amplitude: float, width: float, center: Optional[float] = None
494    ) -> Array:
495        return self.gaussian_transverse_displacement(amplitude, width, center)
496
497
498def cantilever_model(
499    length: float = 4.0,
500    height: float = 1.0,
501    nx: int = 8,
502    ny: int = 2,
503    young_modulus: float = 1.0e6,
504    poisson_ratio: float = 0.3,
505    density: float = 1050.0,
506    mass_type: str = "consistent",
507) -> NonlinearSolid2D:
508    mesh = StructuredQ4Mesh(length, height, nx, ny)
509    material = Material(young_modulus, poisson_ratio, density)
510    constrained = NonlinearSolid2D.clamp_dofs(mesh, left=True, right=False)
511    return NonlinearSolid2D(mesh, material, constrained, mass_type=mass_type)
512
513
514def doubly_clamped_model(
515    length: float = 4.0,
516    height: float = 1.0,
517    nx: int = 16,
518    ny: int = 4,
519    young_modulus: float = 1.0e6,
520    poisson_ratio: float = 0.3,
521    density: float = 1050.0,
522    mass_type: str = "consistent",
523) -> NonlinearSolid2D:
524    mesh = StructuredQ4Mesh(length, height, nx, ny)
525    material = Material(young_modulus, poisson_ratio, density)
526    constrained = NonlinearSolid2D.clamp_dofs(mesh, left=True, right=True)
527    return NonlinearSolid2D(mesh, material, constrained, mass_type=mass_type)

Material-model facade and linear elasticity:

  1"""Material-model facade for the lightweight solid-dynamics examples.
  2
  3The nonlinear implementation lives in :mod:`nonlinear_solid_dynamics`.  This
  4module adds a small-strain isotropic linear-elastic material while reusing the
  5same mesh, mass matrix, boundary-condition handling, and explicit/implicit time
  6integrators.
  7"""
  8
  9from __future__ import annotations
 10
 11from typing import Callable, Optional
 12
 13import numpy as np
 14
 15try:
 16    import nonlinear_solid_dynamics as nonlinear
 17except ImportError:  # pragma: no cover
 18    from . import nonlinear_solid_dynamics as nonlinear
 19
 20Array = np.ndarray
 21Material = nonlinear.Material
 22SolidState = nonlinear.SolidState
 23StructuredQ4Mesh = nonlinear.StructuredQ4Mesh
 24NonlinearSolid2D = nonlinear.NonlinearSolid2D
 25
 26
 27class LinearElasticSolid2D(NonlinearSolid2D):
 28    """Small-strain isotropic linear-elastic Q4 solid in plane strain."""
 29
 30    def _elasticity_tensor(self) -> Array:
 31        lam = self.material.lame_lambda
 32        mu = self.material.lame_mu
 33        eye = np.eye(2)
 34        return (
 35            lam * np.einsum("ij,kl->ijkl", eye, eye)
 36            + mu * np.einsum("ik,jl->ijkl", eye, eye)
 37            + mu * np.einsum("il,jk->ijkl", eye, eye)
 38        )
 39
 40    def _cauchy_stress(self, grad_u: Array) -> Array:
 41        strain = 0.5 * (grad_u + grad_u.T)
 42        return (
 43            self.material.lame_lambda * np.trace(strain) * np.eye(2)
 44            + 2.0 * self.material.lame_mu * strain
 45        )
 46
 47    def _element_force_and_tangent(
 48        self, element_id: int, u_local: Array
 49    ) -> tuple[Array, Array]:
 50        u_local = np.asarray(u_local, dtype=float).reshape(4, 2)
 51        fe = np.zeros((4, 2), dtype=float)
 52        ke4 = np.zeros((4, 2, 4, 2), dtype=float)
 53        elasticity = self._elasticity_tensor()
 54        for _n, grad_n, dv in self._element_quadrature(element_id):
 55            grad_u = u_local.T @ grad_n
 56            stress = self._cauchy_stress(grad_u)
 57            fe += (grad_n @ stress.T) * dv
 58            ke4 += (
 59                np.einsum("aJ,iJkL,bL->aibk", grad_n, elasticity, grad_n) * dv
 60            )
 61        return fe.reshape(-1), ke4.reshape(8, 8)
 62
 63    def element_internal_force(self, element_id: int, u_local: Array) -> Array:
 64        u_local = np.asarray(u_local, dtype=float).reshape(4, 2)
 65        fe = np.zeros((4, 2), dtype=float)
 66        for _n, grad_n, dv in self._element_quadrature(element_id):
 67            grad_u = u_local.T @ grad_n
 68            stress = self._cauchy_stress(grad_u)
 69            fe += (grad_n @ stress.T) * dv
 70        return fe.reshape(-1)
 71
 72    def acceleration(
 73        self,
 74        displacement: Array,
 75        velocity: Array,
 76        time: float,
 77        external_force: Callable[[float], Array],
 78    ) -> Array:
 79        """Evaluate acceleration, using diagonal inversion for lumped mass."""
 80        if self.mass_type != "lumped":
 81            return super().acceleration(displacement, velocity, time, external_force)
 82        c = self.damping_matrix()
 83        rhs = (
 84            external_force(time)
 85            - self.internal_force(displacement)
 86            - c @ velocity
 87        )
 88        acceleration = np.zeros(self.ndof, dtype=float)
 89        mass_diagonal = np.diag(self._mass)
 90        acceleration[self.free_dofs] = (
 91            rhs[self.free_dofs] / mass_diagonal[self.free_dofs]
 92        )
 93        return acceleration
 94
 95
 96def _solid_class(material_model: str):
 97    key = material_model.lower().replace("-", "_")
 98    if key in ("neo_hookean", "neohookean", "nonlinear"):
 99        return NonlinearSolid2D
100    if key in ("linear", "linear_elastic", "linear_elasticity"):
101        return LinearElasticSolid2D
102    raise ValueError(
103        "material_model must be 'neo_hookean' or 'linear' "
104        f"(received {material_model!r})"
105    )
106
107
108def cantilever_model(
109    length: float = 4.0,
110    height: float = 1.0,
111    nx: int = 8,
112    ny: int = 2,
113    young_modulus: float = 1.0e6,
114    poisson_ratio: float = 0.3,
115    density: float = 1050.0,
116    mass_type: str = "consistent",
117    material_model: str = "neo_hookean",
118):
119    """Construct a left-clamped beam with the requested constitutive model."""
120    mesh = StructuredQ4Mesh(length, height, nx, ny)
121    material = Material(young_modulus, poisson_ratio, density)
122    cls = _solid_class(material_model)
123    constrained = cls.clamp_dofs(mesh, left=True, right=False)
124    return cls(mesh, material, constrained, mass_type=mass_type)
125
126
127def doubly_clamped_model(
128    length: float = 4.0,
129    height: float = 1.0,
130    nx: int = 16,
131    ny: int = 4,
132    young_modulus: float = 1.0e6,
133    poisson_ratio: float = 0.3,
134    density: float = 1050.0,
135    mass_type: str = "consistent",
136    material_model: str = "neo_hookean",
137):
138    """Construct a beam with both ends clamped."""
139    mesh = StructuredQ4Mesh(length, height, nx, ny)
140    material = Material(young_modulus, poisson_ratio, density)
141    cls = _solid_class(material_model)
142    constrained = cls.clamp_dofs(mesh, left=True, right=True)
143    return cls(mesh, material, constrained, mass_type=mass_type)
144
145
146def gaussian_displacement(
147    model: NonlinearSolid2D,
148    amplitude: float,
149    width: float,
150    center: Optional[float] = None,
151    direction: str = "transverse",
152) -> Array:
153    """Create an axial or transverse Gaussian nodal displacement."""
154    if center is None:
155        center = 0.5 * model.mesh.length
156    x = model.mesh.coordinates[:, 0]
157    profile = float(amplitude) * np.exp(-0.5 * ((x - center) / float(width)) ** 2)
158    displacement = np.zeros(model.ndof)
159    key = direction.lower()
160    if key in ("axial", "x", "longitudinal"):
161        displacement[0::2] = profile
162    elif key in ("transverse", "y"):
163        displacement[1::2] = profile
164    else:
165        raise ValueError("direction must be 'axial' or 'transverse'")
166    displacement[model.constrained_dofs] = 0.0
167    return displacement
168
169
170def longitudinal_wave_model(
171    length: float = 4.0,
172    height: float = 1.0,
173    nx: int = 80,
174    ny: int = 1,
175    young_modulus: float = 1.0e6,
176    poisson_ratio: float = 0.3,
177    density: float = 1050.0,
178    mass_type: str = "lumped",
179) -> LinearElasticSolid2D:
180    """Construct a 1-D-like plane-strain bar for a scalar wave-equation check.
181
182    All transverse DOFs are constrained and the axial DOFs are fixed at both
183    ends.  For fields that are uniform through the height, the continuum
184    equations reduce to ``u_tt = c_p**2 u_xx`` with
185    ``c_p = sqrt((lambda + 2*mu)/rho)``.
186    """
187    mesh = StructuredQ4Mesh(length, height, nx, ny)
188    material = Material(young_modulus, poisson_ratio, density)
189    transverse = np.arange(1, mesh.ndof, 2, dtype=int)
190    left = mesh.nodes_on_x(0.0)
191    right = mesh.nodes_on_x(mesh.length)
192    axial_ends = 2 * np.concatenate((left, right))
193    constrained = np.unique(np.concatenate((transverse, axial_ends)))
194    return LinearElasticSolid2D(mesh, material, constrained, mass_type=mass_type)
195
196
197def longitudinal_wave_speed(material: Material) -> float:
198    """Plane-strain longitudinal wave speed for ``longitudinal_wave_model``."""
199    return float(
200        np.sqrt((material.lame_lambda + 2.0 * material.lame_mu) / material.density)
201    )
202
203
204def two_way_gaussian_solution(
205    x: Array,
206    time: float,
207    amplitude: float,
208    width: float,
209    center: float,
210    wave_speed: float,
211) -> Array:
212    """D'Alembert solution for a stationary Gaussian displacement pulse.
213
214    This whole-line solution is the exact solution of ``u_tt = c**2 u_xx`` for
215    ``u(x, 0) = g(x)`` and ``u_t(x, 0) = 0``.  It can be compared directly to
216    the clamped-bar solution before the two traveling pulses reach the ends,
217    provided the initial Gaussian is negligible at the boundaries.
218    """
219    x = np.asarray(x, dtype=float)
220    sigma = float(width)
221    c = float(wave_speed)
222    t = float(time)
223    left = np.exp(-0.5 * ((x + c * t - center) / sigma) ** 2)
224    right = np.exp(-0.5 * ((x - c * t - center) / sigma) ** 2)
225    return 0.5 * float(amplitude) * (left + right)

romtools QoI wrapper:

  1"""romtools QoI-model wrapper for the solid-dynamics example."""
  2
  3from __future__ import annotations
  4
  5import os
  6from typing import Dict
  7
  8import numpy as np
  9
 10try:
 11    import solid_dynamics
 12except ImportError:  # pragma: no cover
 13    from . import solid_dynamics
 14
 15
 16class NonlinearSolidBeamQoiModel:
 17    """Parameterised cantilever model with tip displacement as the QoI.
 18
 19    The parameter sample accepts ``young_modulus``, ``load_amplitude``, and
 20    ``pulse_duration``.  The full displacement and velocity snapshot histories
 21    are retained so the same evaluations can later be reused for ROM training.
 22    ``material_model`` may be set to ``"neo_hookean"`` (default) or ``"linear"``
 23    when constructing the wrapper.
 24    """
 25
 26    def __init__(
 27        self,
 28        nx: int = 8,
 29        ny: int = 2,
 30        length: float = 4.0,
 31        height: float = 1.0,
 32        poisson_ratio: float = 0.3,
 33        density: float = 1050.0,
 34        dt: float = 1.0e-2,
 35        t_end: float = 0.6,
 36        snapshot_stride: int = 2,
 37        material_model: str = "neo_hookean",
 38    ) -> None:
 39        self.nx = int(nx)
 40        self.ny = int(ny)
 41        self.length = float(length)
 42        self.height = float(height)
 43        self.poisson_ratio = float(poisson_ratio)
 44        self.density = float(density)
 45        self.dt = float(dt)
 46        self.t_end = float(t_end)
 47        self.snapshot_stride = int(snapshot_stride)
 48        self.material_model = str(material_model)
 49
 50    def populate_run_directory(self, run_directory: str, parameter_sample: Dict) -> None:
 51        os.makedirs(run_directory, exist_ok=True)
 52        with open(os.path.join(run_directory, "params.txt"), "w", encoding="utf-8") as handle:
 53            for key, value in parameter_sample.items():
 54                handle.write(f"{key}: {value}\n")
 55
 56    def run_model(self, run_directory: str, parameter_sample: Dict) -> int:
 57        model = solid_dynamics.cantilever_model(
 58            length=self.length,
 59            height=self.height,
 60            nx=self.nx,
 61            ny=self.ny,
 62            young_modulus=float(parameter_sample["young_modulus"]),
 63            poisson_ratio=self.poisson_ratio,
 64            density=self.density,
 65            mass_type="consistent",
 66            material_model=self.material_model,
 67        )
 68        load_amplitude = float(parameter_sample["load_amplitude"])
 69        pulse_duration = float(parameter_sample["pulse_duration"])
 70        edge_shape = model.boundary_force_x(model.mesh.length, np.array([0.0, -1.0]))
 71
 72        def external_force(time: float) -> np.ndarray:
 73            if 0.0 <= time <= pulse_duration:
 74                scale = load_amplitude * np.sin(np.pi * time / pulse_duration)
 75            else:
 76                scale = 0.0
 77            return scale * edge_shape
 78
 79        state = model.initial_state(external_force=external_force)
 80        num_steps = int(round(self.t_end / self.dt))
 81        times, displacements, velocities = model.solve_implicit(
 82            state,
 83            dt=self.dt,
 84            num_steps=num_steps,
 85            external_force=external_force,
 86            snapshot_stride=self.snapshot_stride,
 87        )
 88        tip_node = model.mesh.node_nearest(self.length, 0.5 * self.height)
 89        qoi = displacements[:, 2 * tip_node + 1]
 90        np.savez(
 91            os.path.join(run_directory, "solution.npz"),
 92            times=times,
 93            displacements=displacements,
 94            velocities=velocities,
 95            qoi=qoi,
 96        )
 97        return 0
 98
 99    def compute_qoi(self, run_directory: str, parameter_sample: Dict) -> np.ndarray:
100        del parameter_sample
101        return np.load(os.path.join(run_directory, "solution.npz"))["qoi"]
102
103
104def _main() -> None:
105    run_directory = os.path.join(os.path.dirname(__file__), "nonlinear_solid_output")
106    model = NonlinearSolidBeamQoiModel()
107    sample = {
108        "young_modulus": 1.0e6,
109        "load_amplitude": 1000.0,
110        "pulse_duration": 0.15,
111    }
112    model.populate_run_directory(run_directory, sample)
113    model.run_model(run_directory, sample)
114    qoi = model.compute_qoi(run_directory, sample)
115    print(f"Saved {os.path.join(run_directory, 'solution.npz')}")
116    print(f"QoI snapshots: {qoi.size}; minimum tip displacement: {qoi.min():.6e} m")
117
118
119if __name__ == "__main__":
120    _main()