[verification] Add verification cases for the LinearSmallStrainFEMForceField

.gitignore

@@ -0,0 +1,5 @@
+# ignore cached bytecode-compiled python scene files
+*__pycache__/
+
+# a run's own output
+verification/results/

README.md

@@ -1,2 +1,45 @@
 # Verification-Validation
 Initially funded by TIRREX, this repository proposes a framework for verification and validation. It relies on both manufactured solutions and comparison with alternative simulation technologies. This work therefore proves the reliability and compliance of the SOFA models & algorithms.
+
+## Verification
+
+**Verification** uses the method of manufactured solutions: an exact displacement is chosen, the
+loads that produce it are derived, and SOFA solves the problem on a sequence of refined meshes.
+The error must shrink at the rate the discretisation predicts. Covered force fields:
+- `LinearSmallStrainFEMForceField`
+
+To get started, read [`verification/docs/IntroductionMMS.ipynb`](verification/docs/IntroductionMMS.ipynb): it presents the method and builds one verification study step by step, from the manufactured solution to the input file the suite runs.
+
+The following force field components are currently verified within this repository:
+- `LinearSmallStrainFEMForceField` [`verification/docs/LinearSmallStrainFEMForceField.ipynb`](verification/docs/LinearSmallStrainFEMForceField.ipynb)
+
+## Layout
+
+```
+common/                                  shared code: geometry and SOFA helpers
+verification/                            the verification suite
+  cases/<Component>/**.json              one test case per file
+```
+
+## Requirements
+
+- SOFA with SofaPython3
+- numpy, sympy, matplotlib
+
+## Running
+
+```bash
+python verification/run.py <path-to-json-file>  # one test case
+python verification/run.py --all                # every test case, then one summary line per case
+python verification/run.py --clean              # delete previous results
+```
+
+Results are written to `verification/results/`.
+
+To inspect one mesh of a case in the SOFA GUI, the runner also acts as a SOFA scene file:
+
+```bash
+runSofa -l SofaPython3 verification/run.py --argv <path-to-json-file> --argv 0
+```
+
+The first argument is the path to the .json case file and the second is the mesh refinement level to be built.

common/geometry/__init__.py

@@ -0,0 +1,3 @@
+from .geometry import Geometry
+from .grid import Grid, Bar1D, Beam2D, Beam3D
+from .tire import Tire, Tire2D, Tire3D

common/geometry/geometry.py

@@ -0,0 +1,116 @@
+"""Tool-agnostic geometry descriptions for the VnV suite."""
+
+from abc import ABC, abstractmethod
+
+
+def _region_indices_from_mesh(mesh_path):
+    """Node indices tagged by each physical group of a Gmsh .msh file, keyed by group name."""
+    import meshio
+
+    mesh = meshio.read(mesh_path)
+    tag_names = {tag: name for name, (tag, _dim) in mesh.field_data.items()}
+    groups = {}
+    for block, tags in zip(mesh.cells, mesh.cell_data.get("gmsh:physical", [])):
+        for element, tag in zip(block.data, tags):
+            name = tag_names.get(tag)
+            if name is not None:
+                groups.setdefault(name, set()).update(element.tolist())
+    return {name: sorted(indices) for name, indices in groups.items()}
+
+
+class Geometry(ABC):
+    """A shape parameterization plus its named boundary regions."""
+
+    dim: int  # topological dimension of the shape (and of the elements meshing it)
+    PARAMETER_ORDER: tuple  # the family's parameters, ordered by the dimension that introduces each
+
+    def __init__(self, *, dim=None, spatialDimensions=None, meshPath=None, **parameters):
+        """Build this shape, taking each name in its `PARAMETER_NAMES` as a keyword.
+
+        e.g. ``Beam3D(length=2.0, width=1.0, height=1.0)``. Pass `dim` only for a family base,
+        which names no dimension of its own: ``Grid(dim=3, length=2.0, width=1.0, height=1.0)``
+        builds that same box.
+        """
+        stated = getattr(type(self), "dim", None)
+        if stated is None and dim is None:
+            raise TypeError(f"{type(self).__name__} states no dim: pass dim= to build one, or "
+                            f"name the geometry that states it")
+        if stated is not None and dim is not None and dim != stated:
+            raise TypeError(f"{type(self).__name__} is {stated}D and cannot be built with dim={dim}")
+        self.dim = stated if dim is None else dim
+
+        if set(parameters) != set(self.PARAMETER_NAMES):
+            raise TypeError(f"{type(self).__name__} takes {list(self.PARAMETER_NAMES)}, "
+                            f"got {sorted(parameters)}")
+
+        for name in self.PARAMETER_NAMES:
+            setattr(self, name, parameters[name])
+        self._set_spatial_dimensions(spatialDimensions)
+
+        if meshPath is not None:
+            self._regions = {name: (lambda nodes, indices=indices: indices)
+                             for name, indices in _region_indices_from_mesh(meshPath).items()}
+        else:
+            self._regions = {name: (lambda nodes, predicate=predicate:
+                                    [i for i, p in enumerate(nodes) if predicate(p)])
+                             for name, predicate in self._analytic_predicates.items()}
+
+    @property
+    def PARAMETER_NAMES(self) -> tuple:
+        """The parameters this shape takes: the first `dim` of its family's.
+
+        Each dimension adds one, so a member takes as many as it has dimensions: of Grid's
+        ("length", "width", "height"), Bar1D takes ("length",) and Beam2D ("length", "width").
+        A family that does not grow this way sets `PARAMETER_NAMES` itself instead.
+        """
+        return self.PARAMETER_ORDER[:self.dim]
+
+    @property
+    def parameters(self) -> list:
+        """This shape's defining parameters, read off the attributes named in `PARAMETER_NAMES`.
+
+        e.g. Beam2D(length=2.0, width=0.5).parameters == [2.0, 0.5]
+        """
+        return [getattr(self, name) for name in self.PARAMETER_NAMES]
+
+    @property
+    def named_parameters(self) -> dict:
+        """Parameters keyed by name: how a manufactured solution asks for one of them.
+
+        e.g. Beam2D(length=2.0, width=0.5).named_parameters == {"length": 2.0, "width": 0.5}
+        """
+        return dict(zip(self.PARAMETER_NAMES, self.parameters))
+
+    def _set_spatial_dimensions(self, spatial_dimensions):
+        """Set how many dimensions of space this shape is embedded in; defaults to its own."""
+        if spatial_dimensions is None:
+            spatial_dimensions = self.dim
+
+        if spatial_dimensions < self.dim:
+            raise ValueError(
+                f"{type(self).__name__} is {self.dim}D and cannot be embedded in a "
+                f"{spatial_dimensions}D space: the embedding space must be at least as "
+                f"large as the shape."
+            )
+        if spatial_dimensions > 3:
+            raise ValueError(
+                f"{type(self).__name__}: spatial_dimensions cannot exceed 3, "
+                f"got {spatial_dimensions}."
+            )
+
+        self.spatial_dimensions = spatial_dimensions
+
+    @property
+    def region_names(self):
+        """This geometry's named boundary regions."""
+        return self._regions.keys()
+
+    def region_indices(self, region, nodes) -> list:
+        """Indices into `nodes` belonging to `region`."""
+        return self._regions[region](nodes)
+
+    @property
+    @abstractmethod
+    def _analytic_predicates(self) -> dict:
+        """Named boundary regions as ``{name: predicate(point) -> bool}``, from this geometry's
+        own shape parameters."""

common/geometry/grid.py

@@ -0,0 +1,36 @@
+"""Rectilinear structured grid geometry in 1 to 3 dimensions: bars and beams."""
+
+import numpy as np
+
+from .geometry import Geometry
+
+
+class Grid(Geometry):
+    """A rectilinear structured grid, axis i spanning 0..parameters[i]. Always starts at the origin."""
+
+    PARAMETER_ORDER = ("length", "width", "height")
+    AXIS_NAMES = (("left", "right"), ("bottom", "top"), ("front", "back"))
+
+    @property
+    def _analytic_predicates(self) -> dict:
+        predicates = {}
+        for axis, (low, high) in enumerate(self.AXIS_NAMES[:self.dim]):
+            extent = self.parameters[axis]
+            predicates[low] = lambda p, axis=axis: np.isclose(p[axis], 0.0)
+            predicates[high] = lambda p, axis=axis, extent=extent: np.isclose(p[axis], extent)
+        return predicates
+
+
+class Bar1D(Grid):
+    """A 1D bar occupying ``[0, length]`` along x."""
+    dim = 1
+
+
+class Beam2D(Grid):
+    """A 2D rectangle occupying ``[0, length] x [0, width]`` in the xy-plane."""
+    dim = 2
+
+
+class Beam3D(Grid):
+    """A 3D box occupying ``[0, length] x [0, width] x [0, height]``."""
+    dim = 3

common/geometry/tire.py

@@ -0,0 +1,32 @@
+"""Tire geometry (Chamberland et al., 2010, Fig. 1)."""
+
+import numpy as np
+
+from .geometry import Geometry
+
+
+class Tire(Geometry):
+    """An annulus in the xy-plane, centered on the origin, extruded along z when it is 3D."""
+
+    PARAMETER_ORDER = ("inner_radius", "outer_radius", "extrusion_length")
+
+    @property
+    def _analytic_predicates(self) -> dict:
+        predicates = {
+            "inner": lambda p: np.isclose(np.hypot(p[0], p[1]), self.inner_radius),
+            "outer": lambda p: np.isclose(np.hypot(p[0], p[1]), self.outer_radius),
+        }
+        if self.dim == 3:
+            predicates["front"] = lambda p: np.isclose(p[2], 0.0)
+            predicates["back"] = lambda p: np.isclose(p[2], self.extrusion_length)
+        return predicates
+
+
+class Tire2D(Tire):
+    """An annulus occupying ``inner_radius <= r <= outer_radius`` of the xy-plane."""
+    dim = 2
+
+
+class Tire3D(Tire):
+    """A Tire2D extruded along z from 0 to extrusion_length."""
+    dim = 3

common/sofa/__init__.py

@@ -0,0 +1 @@
+from .scene import SofaScene

common/sofa/controllers.py

@@ -0,0 +1,47 @@
+"""Post-init SOFA controllers."""
+
+import Sofa.Core
+
+
+class RegionClamp(Sofa.Core.Controller):
+    """Builds one PartialFixedProjectiveConstraint per boundary region, then applies Dirichlet
+    boundary conditions once the simulation has initialized.
+
+    For each region, the nodes belonging to it are found and handed to its constraint as the indices
+    it fixes. If a prescribed displacement is given, those nodes are also moved from rest to rest +
+    u in the constraint's fixed directions before the constraint clamps them there. rest_position
+    itself is never changed.
+    """
+
+    def __init__(self, geometry, dofs, node, prescribe_displacement_on, vec_type,
+                 displacement=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.geometry = geometry
+        self.dofs = dofs
+        self.displacement = displacement
+
+        # Fix dofs partially as determined by MMS `prescribe_displacement_on`
+        self.region_constraints = []
+        for region, fixed_directions in prescribe_displacement_on.items():
+            # Pad the array with zeroes for embedded dimensions: e.g. [1, 0] -> [1, 0, 0]
+            padded = tuple(fixed_directions) + (0,) * (geometry.spatial_dimensions - len(fixed_directions))
+            constraint = node.addObject('PartialFixedProjectiveConstraint', name='clamp' + region,
+                                  template=vec_type, fixedDirections=list(padded))
+            self.region_constraints.append((region, padded, constraint))
+
+    def onSimulationInitDoneEvent(self, event):
+        rest_positions = self.dofs.rest_position.array()
+
+        # Apply non-zero Dirichlet BC: Move each fixed direction to rest + u; rest_position is untouched.
+        with self.dofs.position.writeableArray() as pos:
+            for region, fixed_directions, constraint in self.region_constraints:
+                indices = self.geometry.region_indices(region, rest_positions)
+
+                if self.displacement is not None and indices:
+                    displacement = self.displacement(rest_positions[indices])
+
+                    for axis, fixed in enumerate(fixed_directions):
+                        if fixed:
+                            pos[indices, axis] = rest_positions[indices, axis] + displacement[:, axis]
+
+                constraint.indices.value = indices

common/sofa/conventions.py

@@ -0,0 +1,38 @@
+"""Central map from tool-agnostic vocabulary to SOFA components."""
+
+from dataclasses import dataclass
+
+# Spatial dimension (DOF/embedding space) -> MechanicalObject / force-field vector template.
+VEC_DIM = {1: "Vec1d", 2: "Vec2d", 3: "Vec3d"}
+
+@dataclass(frozen=True)
+class ElementKind:
+    """Everything the scene-building layer needs to know about one type of finite element."""
+
+    dim: int                    # topological dimension of this element
+    container: str              # topology container class
+    data_name: str              # Data field this kind's own container addresses its elements by
+    cpp: str                    # SOFA geometry name, for compound templates ("Vec3d,Hexahedron")
+    facet_kind: str | None      # the kind of this element's own facets, None if it has none
+    grid_mapping: str | None    # topological mapping generating it from the grid, None if grid-native
+    boundary_mappings: tuple    # (element kind, mapping) steps from this element to its boundary
+
+ELEMENTS = {
+    "edge": ElementKind(dim=1, container="EdgeSetTopologyContainer",        data_name="edges",
+                        cpp="Edge",        facet_kind=None,   grid_mapping=None,
+                        boundary_mappings=()),
+    "tri":  ElementKind(dim=2, container="TriangleSetTopologyContainer",    data_name="triangles",
+                        cpp="Triangle",    facet_kind="edge", grid_mapping="Quad2TriangleTopologicalMapping",
+                        boundary_mappings=(("edge", "Triangle2EdgeTopologicalMapping"),)),
+    "quad": ElementKind(dim=2, container="QuadSetTopologyContainer",        data_name="quads",
+                        cpp="Quad",        facet_kind="edge", grid_mapping=None,
+                        # To arrive at a Quad2Edge map -> triangulate quads and then Tri2Edge
+                        boundary_mappings=(("tri", "Quad2TriangleTopologicalMapping"),
+                                           ("edge", "Triangle2EdgeTopologicalMapping"))),
+    "tet":  ElementKind(dim=3, container="TetrahedronSetTopologyContainer", data_name="tetrahedra",
+                        cpp="Tetrahedron", facet_kind="tri",  grid_mapping="Hexa2TetraTopologicalMapping",
+                        boundary_mappings=(("tri", "Tetra2TriangleTopologicalMapping"),)),
+    "hexa": ElementKind(dim=3, container="HexahedronSetTopologyContainer",  data_name="hexahedra",
+                        cpp="Hexahedron",  facet_kind="quad", grid_mapping=None,
+                        boundary_mappings=(("quad", "Hexa2QuadTopologicalMapping"),)),
+}

common/sofa/prefabs/__init__.py

@@ -0,0 +1,8 @@
+from .base import ScenePrefab
+from .grid import ElasticGrid, ElasticBar, ElasticBeam
+from .tire import ElasticTire
+from ...geometry import Grid, Bar1D, Beam2D, Beam3D, Tire, Tire2D, Tire3D
+
+PREFAB_BY_GEOMETRY = {Grid: ElasticGrid, Bar1D: ElasticBar,
+                      Beam2D: ElasticBeam, Beam3D: ElasticBeam,
+                      Tire: ElasticTire, Tire2D: ElasticTire, Tire3D: ElasticTire}