Adds a complete (MMS) verification framework for 3D tetrahedral (P1) elements

examples/Freefem/MMS/3D/manufactured_solution.py

@@ -20,6 +20,7 @@ class MMSCase3D(MMSCase):
     # Body-force quadrature rule — must be set by each concrete case.
     # No framework fallback; assembly raises if unset.
     source_quadrature_hex = None   # element rule for Q1 hexes (e.g. hex_q1_rule(2))
+    source_quadrature_tet = None   # element rule for P1 tets  (e.g. tet_p1_rule(4))
 
     @abstractmethod
     def u_ex(self, x, y, z, L):

examples/Freefem/MMS/3D/params.json

@@ -3,19 +3,26 @@
     "youngModulus": 1e6,
     "forceField": "LinearSmallStrainFEMForceField",
     "linearSolver": {
-        "type": "SparseLDLSolver",
+        "type": "PCGLinearSolver",
         "parameters": {
-            "template": "CompressedRowSparseMatrixd"
-        }
+            "template": "GraphScattered",
+            "tolerance": 1e-14,
+            "iterations": 200,
+            "preconditioner": "@precond"
+        },
+        "preconditioner": "SSORPreconditioner"
     },
     "reference": {
-        "nx": 6,
+        "nx": 20,
         "nu": 0.3
     },
     "convergence": {
-        "nu_values": [0.3, 0.49],
+        "nu_values": [0.3],
         "nx_values": {
-            "sinusoidal": [10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]
+            "sinusoidal": [
+                  8, 12, 16, 20, 24, 28,
+                32, 36, 40, 44, 48, 52, 56, 60
+            ]
         }
     }
 }

examples/Freefem/MMS/3D/run_convergence.py

@@ -12,6 +12,7 @@
     load_params,
     solve_solid,
     element_hex,
+    element_tet,
 )
 from convergence import run_convergence_series
 from output      import plot_convergence
@@ -70,8 +71,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values,
     conv = cfg["convergence"]
 
     specs = [
-        {"elem": element_hex, "label": "Q1 hex",
-         "l2_style": "bo-", "h1_style": "rs--"},
+        {"elem": element_hex, "label": "Q1 hex", "l2_style": "bo-", "h1_style": "rs--"},
+        {"elem": element_tet, "label": "P1 tet", "l2_style": "g^-", "h1_style": "md--"},
     ]
 
     for mms in (sinusoidal_mms,):

examples/Freefem/MMS/3D/sinusoidal.py

@@ -19,7 +19,8 @@
 
 from manufactured_solution import MMSCase3D, lame
 from solid import (case_scene, run_reference_scene,
-                   element_hex, hex_q1_rule)
+                   element_hex, hex_q1_rule,
+                   element_tet, tet_p1_rule)
 
 
 SINUS_AMPLITUDE = 1e-1
@@ -32,6 +33,7 @@ class SinusNeumann(MMSCase3D):
                   r"\sin(\pi z/L)\sin(\pi x/L))$")
 
     source_quadrature_hex = staticmethod(hex_q1_rule(2))
+    source_quadrature_tet = staticmethod(tet_p1_rule(4))
 
     def u_ex(self, x, y, z, L):
         k = np.pi / L
@@ -84,30 +86,23 @@ def source(self, x, y, z, E, nu, L):
     def apply_bcs(self, Solid, nodes_3d, L):
         eps = 1e-10
         xyz = nodes_3d[:, :3]
+        dofs = Solid.dofs
 
-        def find_corner(pred):
-            for k, (xk, yk, zk) in enumerate(xyz):
-                if pred(xk, yk, zk):
-                    return k
-            raise RuntimeError("sinusoidal: BC corner not found")
+        face = [i for i, (x, y, z) in enumerate(xyz) if x < eps]
 
-        i_origin = find_corner(lambda x, y, z: x < eps     and y < eps     and z < eps)
-        i_xL     = find_corner(lambda x, y, z: x > L - eps and y < eps     and z < eps)
-        i_yL     = find_corner(lambda x, y, z: x < eps     and y > L - eps and z < eps)
+        with dofs.position.writeable() as pos:
+            for i in face:
+                x, y, z = xyz[i]
+                ux, uy, uz = self.u_ex(x, y, z, L)
+                pos[i] = [x + ux, y + uy, z + uz]
 
         Solid.addObject("FixedProjectiveConstraint",
-                        name="fix_origin", indices=i_origin)
-        Solid.addObject("PartialFixedProjectiveConstraint",
-                        name="fix_x_axis", template="Vec3d",
-                        indices=i_xL, fixedDirections="0 1 1")
-        Solid.addObject("PartialFixedProjectiveConstraint",
-                        name="fix_y_axis", template="Vec3d",
-                        indices=i_yL, fixedDirections="0 0 1")
+                        name="fix_face", indices=face)
 
 
 mms         = SinusNeumann()
-createScene = case_scene(mms, element_hex)
+createScene = case_scene(mms, element_tet)
 
 
 if __name__ == "__main__":
-    run_reference_scene(element_hex, mms)
+    run_reference_scene(element_tet, mms)

examples/Freefem/MMS/3D/solid.py

@@ -11,8 +11,8 @@
 
 # Make the parent MMS/ directory importable so we can pull in fem.py.
 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-from fem import hex_q1_rule              # re-exported for case files
-from elements import element_hex          # re-exported for case files
+from fem import hex_q1_rule, tet_p1_rule  # re-exported for case files
+from elements import element_hex, element_tet  # re-exported for case files
 from solid_solution import SolidSolution3D
 from output import write_solution_table
 from scene import NodalForceAssembler
@@ -73,6 +73,8 @@ def build_solid_scene(rootNode, mms, element, force_field, linear_solver,
         "Sofa.Component.Engine.Select",
         "Sofa.Component.LinearSolver.Direct",
         "Sofa.Component.LinearSolver.Iterative",
+        "Sofa.Component.LinearSolver.Preconditioner",
+        "Sofa.Component.LinearSystem",
         "Sofa.Component.MechanicalLoad",
         "Sofa.Component.ODESolver.Backward",
         "Sofa.Component.StateContainer",
@@ -99,20 +101,42 @@ def build_solid_scene(rootNode, mms, element, force_field, linear_solver,
     Solid = rootNode.addChild("Solid")
     Solid.addObject("StaticSolver", name="staticSolver", printLog=False)
     Solid.addObject("NewtonRaphsonSolver", name="newtonSolver",
-                    maxNbIterationsNewton=1,
+                    maxNbIterationsNewton=10,
                     absoluteResidualStoppingThreshold=1e-10,
-                    printLog=False)
-    Solid.addObject(linear_solver["type"], name="linearSolver",
-                    **linear_solver["parameters"])
+                    relativeSuccessiveStoppingThreshold=1e-12,
+                    relativeInitialStoppingThreshold=1e-12,
+                    printLog=True)
+
+    if linear_solver["preconditioner"]:
+        preconditioner = linear_solver["preconditioner"]
+        Solid.addObject("PreconditionedMatrixFreeSystem"
+            , name="solverSystem"
+            , template="GraphScattered"
+            , preconditionerSystem="@precondSystem")
+        Solid.addObject(linear_solver["type"]
+            , name="linearSolver"
+            , **linear_solver["parameters"]
+            , printLog=True)
+
+        Solid.addObject("MatrixLinearSystem"
+            , name="precondSystem"
+            , template="CompressedRowSparseMatrixd")
+        Solid.addObject(preconditioner, name="precond", printLog=True)
+
+    else:
+        Solid.addObject(linear_solver["type"], name="linearSolver", **linear_solver["parameters"])
 
     dofs = Solid.addObject("MechanicalObject", name="dofs", template="Vec3d",
                            position=nodes_3d.tolist(),
+                           rest_position=nodes_3d.tolist(),
                            showObject=with_visual, showObjectScale=0.005 * L)
 
     topology = element.add_topology(Solid)
 
     Solid.addObject(force_field, name="FEM", template="Vec3d",
-                    youngModulus=E, poissonRatio=nu, topology="@topology")
+                    youngModulus=E, poissonRatio=nu, topology="@topology",
+                    computeForceStrategy="parallel", computeForceDerivStrategy="parallel")
+
 
     mms.apply_bcs(Solid, nodes_3d, L)
 
@@ -143,16 +167,15 @@ def solve_solid(elem, mms, L, E, nu, nx, ny, nz, force_field, linear_solver):
         nx=nx, ny=ny, nz=nz, with_visual=False,
         force_field=force_field, linear_solver=linear_solver
     )
-    Sofa.Simulation.init(root)
+    Sofa.Simulation.initRoot(root)
     nodes_3d = dofs.rest_position.array().copy()
     conn     = elem.read_connectivity(topology)
-    pos0     = dofs.position.array().copy()
     Sofa.Simulation.animate(root, root.dt.value)
     pos1     = dofs.position.array().copy()
     Sofa.Simulation.unload(root)
-    ux = pos1[:, 0] - pos0[:, 0]
-    uy = pos1[:, 1] - pos0[:, 1]
-    uz = pos1[:, 2] - pos0[:, 2]
+    ux = pos1[:, 0] - nodes_3d[:, 0]
+    uy = pos1[:, 1] - nodes_3d[:, 1]
+    uz = pos1[:, 2] - nodes_3d[:, 2]
     return SolidSolution3D(nodes=nodes_3d, conn=conn, ux=ux, uy=uy, uz=uz)
 
 

examples/Freefem/MMS/elements.py

@@ -13,10 +13,12 @@
     compute_h1(sol, mms, L)     : H¹ semi-norm error on a solution
 
 Two `_ElementBase` classes live here — one for 2D (4 boundary edges,
-dim-dependent plane stress / plane strain) and one for 3D (6 boundary
-quads, single constitutive branch). Concrete elements (`_QuadElement`,
-`_TriElement`, `_HexElement`) subclass the matching base and pin the
-topology container choice and the rules.
+dim-dependent plane stress / plane strain) and one for 3D (six domain
+faces, single constitutive branch). The 3D base delegates the boundary
+facets to `_boundary_facet_groups` / `_facet_rule` hooks so hex uses grid
+quads and tet uses the triangulated split. Concrete elements
+(`_QuadElement`, `_TriElement`, `_HexElement`, `_TetElement`) subclass the
+matching base and pin the topology container choice and the rules.
 """
 
 import numpy as np
@@ -29,8 +31,10 @@
     quad_q1_rule,
     tri_p1_rule,
     hex_q1_rule,
+    tet_p1_rule,
     edge_line_rule,
     quad_face_rule,
+    tri_face_rule,
 )
 
 
@@ -74,6 +78,59 @@ def idx(i, j, k):
     return xm, xp, ym, yp, zm, zp
 
 
+# Local vertex triples for the four triangular faces of a tetrahedron.
+_TET_LOCAL_FACES = ((1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2))
+
+
+def _boundary_tris_from_conn(conn, nodes):
+    """Boundary triangles of a tet mesh, grouped by the six domain faces.
+    """
+    conn = np.asarray(conn)
+
+    # Count how many tets each triangular face (as a sorted node triple) touches;
+    # keep one representative ordering per face for later coordinate lookup.
+    face_count = {}
+    face_nodes = {}
+    for tet in conn:
+        for a, b, c in _TET_LOCAL_FACES:
+            tri = (int(tet[a]), int(tet[b]), int(tet[c]))
+            key = tuple(sorted(tri))
+            face_count[key] = face_count.get(key, 0) + 1
+            face_nodes.setdefault(key, tri)
+
+    boundary = [face_nodes[k] for k, n in face_count.items() if n == 1]
+
+    xyz = np.asarray(nodes)[:, :3]
+    lo, hi = xyz.min(axis=0), xyz.max(axis=0)
+    tol = 1e-7 * float(np.max(hi - lo))
+
+    planes = {(-1.0, 0.0, 0.0): [], (+1.0, 0.0, 0.0): [],
+              (0.0, -1.0, 0.0): [], (0.0, +1.0, 0.0): [],
+              (0.0, 0.0, -1.0): [], (0.0, 0.0, +1.0): []}
+    for tri in boundary:
+        p = xyz[list(tri)]                        # (3, 3)
+        matched = False
+        for axis in range(3):
+            if np.all(np.abs(p[:, axis] - lo[axis]) < tol):
+                nrm = [0.0, 0.0, 0.0]; nrm[axis] = -1.0
+            elif np.all(np.abs(p[:, axis] - hi[axis]) < tol):
+                nrm = [0.0, 0.0, 0.0]; nrm[axis] = +1.0
+            else:
+                continue
+            planes[tuple(nrm)].append(tri)
+            matched = True
+            break
+        if not matched:
+            # A boundary face that lies on no domain plane means the mesh is
+            # not the expected [0,L]^3 box tessellation
+            raise RuntimeError(
+                "boundary triangle not on any domain face plane; "
+                "unexpected mesh geometry")
+
+    return [(np.asarray(tris, dtype=int), nrm)
+            for nrm, tris in planes.items() if tris]
+
+
 # ---------------------------------------------------------------------------
 # 2D elements (Vec2d / plane stress  OR  Vec3d / plane strain via dim arg)
 # ---------------------------------------------------------------------------
@@ -192,19 +249,13 @@ def compute_nodal_forces(cls, nodes_3d, conn, mms, L, E, nu, nx, ny, nz):
             lambda x, y, z: mms.source(x, y, z, E, nu, L),
             xyz, conn, cls._source_rule(mms))
 
-        xm, xp, ym, yp, zm, zp = _boundary_quads(nx, ny, nz)
-        sides = [(xm, -1.0, 0.0, 0.0),
-                 (xp, +1.0, 0.0, 0.0),
-                 (ym,  0.0, -1.0, 0.0),
-                 (yp,  0.0, +1.0, 0.0),
-                 (zm,  0.0, 0.0, -1.0),
-                 (zp,  0.0, 0.0, +1.0)]
-        face_rule = quad_face_rule(2)
-        for quads, nrm_x, nrm_y, nrm_z in sides:
+        facet_rule = cls._facet_rule()
+        for facets, (nrm_x, nrm_y, nrm_z) in cls._boundary_facet_groups(
+                conn, xyz, nx, ny, nz):
             F += assemble_traction(
                 lambda x, y, z, nx=nrm_x, ny=nrm_y, nz=nrm_z:
                     mms.traction(x, y, z, nx, ny, nz, E, nu, L),
-                xyz, quads, face_rule)
+                xyz, facets, facet_rule)
         return F
 
     @classmethod
@@ -236,6 +287,20 @@ def _source_rule(mms):
                 f"{type(mms).__name__}.source_quadrature_hex must be set")
         return rule
 
+    @staticmethod
+    def _facet_rule():
+        return quad_face_rule(2)
+
+    @staticmethod
+    def _boundary_facet_groups(conn, xyz, nx, ny, nz):
+        xm, xp, ym, yp, zm, zp = _boundary_quads(nx, ny, nz)
+        return [(xm, (-1.0, 0.0, 0.0)),
+                (xp, (+1.0, 0.0, 0.0)),
+                (ym, (0.0, -1.0, 0.0)),
+                (yp, (0.0, +1.0, 0.0)),
+                (zm, (0.0, 0.0, -1.0)),
+                (zp, (0.0, 0.0, +1.0))]
+
     @staticmethod
     def add_topology(Solid):
         topology = Solid.addObject("HexahedronSetTopologyContainer",
@@ -250,10 +315,48 @@ def read_connectivity(topology):
         return topology.hexahedra.array().copy()
 
 
+class _TetElement(_ElementBase3D):
+    LABEL        = "P1 tet"
+    ELEMENT_RULE = staticmethod(tet_p1_rule(4))   # used for L²/H¹ error norms
+
+    @staticmethod
+    def _source_rule(mms):
+        rule = mms.source_quadrature_tet
+        if rule is None:
+            raise ValueError(
+                f"{type(mms).__name__}.source_quadrature_tet must be set")
+        return rule
+
+    @staticmethod
+    def _facet_rule():
+        return tri_face_rule(3)
+
+    @staticmethod
+    def _boundary_facet_groups(conn, xyz, nx, ny, nz):
+        # Tet boundary faces are the true once-used triangles of the actual
+        # split — independent of the Hexa2Tetra diagonal orientation.
+        return _boundary_tris_from_conn(conn, xyz)
+
+    @staticmethod
+    def add_topology(Solid):
+        topology = Solid.addObject("TetrahedronSetTopologyContainer",
+                                   name="topology")
+        Solid.addObject("Hexa2TetraTopologicalMapping",
+                        input="@../Grid/grid", output="@topology",
+                        swapping=True)
+        Solid.addObject("TetrahedronSetTopologyModifier")
+        return topology
+
+    @staticmethod
+    def read_connectivity(topology):
+        return topology.tetrahedra.array().copy()
+
+
 # ---------------------------------------------------------------------------
 # Instances (one per element type)
 # ---------------------------------------------------------------------------
 
 element_quad = _QuadElement()
 element_tri  = _TriElement()
 element_hex  = _HexElement()
+element_tet  = _TetElement()