diff --git a/examples/Freefem/MMS/3D/manufactured_solution.py b/examples/Freefem/MMS/3D/manufactured_solution.py index 6d64839f..fea38ad9 100644 --- a/examples/Freefem/MMS/3D/manufactured_solution.py +++ b/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): diff --git a/examples/Freefem/MMS/3D/params.json b/examples/Freefem/MMS/3D/params.json index 3e3652e6..bacf2faa 100644 --- a/examples/Freefem/MMS/3D/params.json +++ b/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 + ] } } } diff --git a/examples/Freefem/MMS/3D/run_convergence.py b/examples/Freefem/MMS/3D/run_convergence.py index 16a1e6ae..eb86dc7a 100644 --- a/examples/Freefem/MMS/3D/run_convergence.py +++ b/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,): diff --git a/examples/Freefem/MMS/3D/sinusoidal.py b/examples/Freefem/MMS/3D/sinusoidal.py index f057552e..7188b84f 100644 --- a/examples/Freefem/MMS/3D/sinusoidal.py +++ b/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) diff --git a/examples/Freefem/MMS/3D/solid.py b/examples/Freefem/MMS/3D/solid.py index 03c01fcb..a6050f02 100644 --- a/examples/Freefem/MMS/3D/solid.py +++ b/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) diff --git a/examples/Freefem/MMS/elements.py b/examples/Freefem/MMS/elements.py index 1f4fc9c4..cb6d8781 100644 --- a/examples/Freefem/MMS/elements.py +++ b/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,6 +315,43 @@ 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) # --------------------------------------------------------------------------- @@ -257,3 +359,4 @@ def read_connectivity(topology): element_quad = _QuadElement() element_tri = _TriElement() element_hex = _HexElement() +element_tet = _TetElement() diff --git a/examples/Freefem/MMS/fem.py b/examples/Freefem/MMS/fem.py index edb7e0b1..7aae6f60 100644 --- a/examples/Freefem/MMS/fem.py +++ b/examples/Freefem/MMS/fem.py @@ -7,8 +7,8 @@ `h1_semi_error_1d`. * **Dim-agnostic vector path** used by the 2D and 3D drivers: - `quad_q1_rule`, `tri_p1_rule`, `hex_q1_rule` (element rules), - `edge_line_rule`, `quad_face_rule` (facet rules), + `quad_q1_rule`, `tri_p1_rule`, `hex_q1_rule`, `tet_p1_rule` (element rules), + `edge_line_rule`, `quad_face_rule`, `tri_face_rule` (facet rules), `assemble_nodal_forces`, `assemble_traction`, `l2_error`, `h1_semi_error`. @@ -129,6 +129,22 @@ def _shape_hex_q1(xi, eta, zeta): } +# Reference-tetrahedron quadrature over {(xi,eta,zeta): >=0, xi+eta+zeta<=1}. +# Weights sum to the reference volume 1/6; the physical weight is w_ref * 6*vol. +# The 4-point rule is degree-2 exact (a = (5-sqrt5)/20, b = (5+3sqrt5)/20). +_TET_A4 = (5.0 - np.sqrt(5.0)) / 20.0 +_TET_B4 = (5.0 + 3.0 * np.sqrt(5.0)) / 20.0 +_TET_QUADRATURE = { + 1: (np.array([[1/4, 1/4, 1/4]]), + np.array([1/6])), + 4: (np.array([[_TET_A4, _TET_A4, _TET_A4], + [_TET_B4, _TET_A4, _TET_A4], + [_TET_A4, _TET_B4, _TET_A4], + [_TET_A4, _TET_A4, _TET_B4]]), + np.array([1/24, 1/24, 1/24, 1/24])), +} + + # --------------------------------------------------------------------------- # Element rules (consume by assemble_nodal_forces, l2_error, h1_semi_error) # --------------------------------------------------------------------------- @@ -211,6 +227,35 @@ def rule(xe): return rule +def tet_p1_rule(n_pts=4): + """Element rule for P1 tetrahedra: 1-point (centroid) or 4-point Gauss. + """ + if n_pts not in _TET_QUADRATURE: + raise ValueError(f"tet_p1_rule: {n_pts}-point rule not supported") + pts, wts = _TET_QUADRATURE[n_pts] + + # Reference-coord gradients of (N0..N3) wrt (xi, eta, zeta), constant. + dN_ref = np.array([[-1.0, -1.0, -1.0], + [ 1.0, 0.0, 0.0], + [ 0.0, 1.0, 0.0], + [ 0.0, 0.0, 1.0]]) # (4, 3) [a, ref-axis] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (4, 3) + x0 = xe_arr[0] + # Jacobian columns are the edge vectors from node 0: x = x0 + J @ (xi,eta,zeta). + J = np.column_stack([xe_arr[1] - x0, xe_arr[2] - x0, xe_arr[3] - x0]) + detJ = np.linalg.det(J) + vol = abs(detJ) / 6.0 + Jinv = np.linalg.inv(J) + dN_phys = (dN_ref @ Jinv).T # (3, 4) [d, a] + for (xi, eta, zeta), w_ref in zip(pts, wts): + N = np.array([1.0 - xi - eta - zeta, xi, eta, zeta]) + coords = N @ xe_arr # (3,) + yield coords, w_ref * 6.0 * vol, N, dN_phys + return rule + + # --------------------------------------------------------------------------- # Facet rules (consumed by assemble_traction) # --------------------------------------------------------------------------- @@ -265,6 +310,27 @@ def rule(xe): return rule +def tri_face_rule(n_pts=3): + """Facet rule for a 3D boundary triangle: linear shape, reference-triangle Gauss. + + Yields per Gauss point on the face: (coords (3,), w, N (3,)). + Area via dA = ||(x1-x0) x (x2-x0)|| / 2; vertex order does not affect it. + """ + if n_pts not in _TRI_QUADRATURE: + raise ValueError(f"tri_face_rule: {n_pts}-point rule not supported") + pts, wts = _TRI_QUADRATURE[n_pts] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (3, 3) + x0, x1, x2 = xe_arr + area = 0.5 * np.linalg.norm(np.cross(x1 - x0, x2 - x0)) + for (xi, eta), w_ref in zip(pts, wts): + N = np.array([1.0 - xi - eta, xi, eta]) + coords = N @ xe_arr # (3,) + yield coords, w_ref * 2.0 * area, N + return rule + + # --------------------------------------------------------------------------- # 1D FEM assembly (scalar-integrand path) # ---------------------------------------------------------------------------