[3D] Rethink 3D problem with tetrahedra

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)