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/3D_tetra/convergence.py b/examples/Freefem/MMS/3D_tetra/convergence.py deleted file mode 100644 index ddf98824..00000000 --- a/examples/Freefem/MMS/3D_tetra/convergence.py +++ /dev/null @@ -1,62 +0,0 @@ -"""Shared convergence-study helpers for MMS 1D/2D/3D drivers.""" - -import numpy as np - -from output import write_convergence_table - - -def run_convergence_series(nx_values, run_fn, h_fn, error_fns, banner, - results_dir, table_stem): - """ - Run a mesh-refinement series and write a table. - - For each `nx` in `nx_values`, runs `run_fn(nx)` to get a solution, - evaluates every error in `error_fns`, prints a live row, and writes - the assembled table to ``/.txt``. - - Parameters - ---------- - nx_values : iterable of mesh-refinement parameters - run_fn : callable(nx) -> Solution (anything `error_fns` understands) - h_fn : callable(nx) -> mesh size h (drives the rate denominator - and the printed h column) - error_fns : ordered dict { label: callable(sol) -> float } - banner : string printed before the table header - results_dir, table_stem : the file ``/.txt`` - receives the assembled table - - Returns - ------- - (hs, errors) : list of mesh sizes; dict {label: list of errors}. - """ - hs = [] - errors = {label: [] for label in error_fns} - rows = [] - err_labels = list(error_fns.keys()) - - hdr = f"{'nx':>5} | {'h':>10}" + "".join( - f" | {lbl:>14} | {('rate_' + lbl):>7}" for lbl in err_labels) - print(f"\n{banner}\n{hdr}", flush=True) - - for k, nx in enumerate(nx_values): - h = h_fn(nx) - sol = run_fn(nx) - row = {"nx": nx, "h": h} - for label, err_fn in error_fns.items(): - e_k = err_fn(sol) - rate = (f"{np.log(e_k / errors[label][-1]) / np.log(h / hs[-1]):.2f}" - if k > 0 else "") - errors[label].append(e_k) - row[label] = e_k - row[f"rate_{label}"] = rate - - line = f"{nx:5d} | {h:10.4f}" - for lbl in err_labels: - line += f" | {row[lbl]:14.6e} | {row[f'rate_{lbl}']:>7}" - print(line, flush=True) - - hs.append(h) - rows.append(row) - - write_convergence_table(table_stem, rows, results_dir) - return hs, errors diff --git a/examples/Freefem/MMS/3D_tetra/elements.py b/examples/Freefem/MMS/3D_tetra/elements.py deleted file mode 100644 index d7f7f75b..00000000 --- a/examples/Freefem/MMS/3D_tetra/elements.py +++ /dev/null @@ -1,351 +0,0 @@ -"""Element strategies for MMS scenes. -... -""" - -import numpy as np - -from fem import ( - assemble_nodal_forces, - assemble_traction, - l2_error, - h1_semi_error, - quad_q1_rule, - tri_p1_rule, - hex_q1_rule, - tet_p1_rule, - edge_line_rule, - quad_face_rule, - tri_face_rule, # <-- AJOUT : nécessaire pour _TetElement -) - - -# --------------------------------------------------------------------------- -# Boundary-facet helpers -# --------------------------------------------------------------------------- - -def _boundary_edges(nx, ny): - bottom = [(i, i + 1) for i in range(nx - 1)] - top = [((ny - 1) * nx + i, (ny - 1) * nx + i + 1) for i in range(nx - 1)] - left = [(j * nx, (j + 1) * nx) for j in range(ny - 1)] - right = [(j * nx + (nx - 1), (j + 1) * nx + (nx - 1)) for j in range(ny - 1)] - return bottom, top, left, right - - -def _boundary_quads(nx, ny, nz): - def idx(i, j, k): - return i + j * nx + k * nx * ny - - xm = [(idx(0, j, k), idx(0, j+1, k), idx(0, j+1, k+1), idx(0, j, k+1)) - for k in range(nz - 1) for j in range(ny - 1)] - xp = [(idx(nx-1, j, k), idx(nx-1, j+1, k), idx(nx-1, j+1, k+1), idx(nx-1, j, k+1)) - for k in range(nz - 1) for j in range(ny - 1)] - ym = [(idx(i, 0, k), idx(i+1, 0, k), idx(i+1, 0, k+1), idx(i, 0, k+1)) - for k in range(nz - 1) for i in range(nx - 1)] - yp = [(idx(i, ny-1, k), idx(i+1, ny-1, k), idx(i+1, ny-1, k+1), idx(i, ny-1, k+1)) - for k in range(nz - 1) for i in range(nx - 1)] - zm = [(idx(i, j, 0), idx(i+1, j, 0), idx(i+1, j+1, 0), idx(i, j+1, 0)) - for j in range(ny - 1) for i in range(nx - 1)] - zp = [(idx(i, j, nz-1), idx(i+1, j, nz-1), idx(i+1, j+1, nz-1), idx(i, j+1, nz-1)) - for j in range(ny - 1) for i in range(nx - 1)] - return xm, xp, ym, yp, zm, zp - - -def _boundary_triangles(nx, ny, nz, diagonal="main"): - """Découpe chaque quad de _boundary_quads en 2 triangles. - - diagonal="main" coupe (0,1,2) & (0,2,3) ; "anti" coupe (1,2,3) & (1,3,0). - Le sens qui matche la vraie diagonale interne de - Hexa2TetraTopologicalMapping(swapping=True) n'est pas garanti a priori - -> il faut tester les deux et garder celui qui stabilise le taux H1. - """ - xm, xp, ym, yp, zm, zp = _boundary_quads(nx, ny, nz) - - def split(quads): - tris = [] - for q in quads: - if diagonal == "main": - tris.append((q[0], q[1], q[2])) - tris.append((q[0], q[2], q[3])) - else: - tris.append((q[1], q[2], q[3])) - tris.append((q[1], q[3], q[0])) - return tris - - return [split(f) for f in (xm, xp, ym, yp, zm, zp)] - - -# --------------------------------------------------------------------------- -# 2D elements (Vec2d / plane stress OR Vec3d / plane strain via dim arg) -# --------------------------------------------------------------------------- - -class _ElementBase2D: - @classmethod - def compute_nodal_forces(cls, nodes_2d, conn, mms, L, E, nu, nx, ny, dim): - xy = nodes_2d[:, :2] - - F = assemble_nodal_forces( - lambda x, y: mms.source(x, y, E, nu, L, dim), - xy, conn, cls._source_rule(mms)) - - bottom, top, left, right = _boundary_edges(nx, ny) - sides = [(bottom, 0.0, -1.0), - (top, 0.0, +1.0), - (left, -1.0, 0.0), - (right, +1.0, 0.0)] - edge_rule = edge_line_rule(2) - for edges, nrm_x, nrm_y in sides: - F += assemble_traction( - lambda x, y, nx=nrm_x, ny=nrm_y: - mms.traction(x, y, nx, ny, E, nu, L, dim), - xy, edges, edge_rule) - return F - - @classmethod - def compute_l2(cls, sol, mms, L): - xy = sol.nodes[:, :2] - return l2_error( - xy, sol.conn, np.column_stack([sol.ux, sol.uy]), - lambda x, y: mms.u_ex(x, y, L), - cls.ELEMENT_RULE) - - @classmethod - def compute_h1(cls, sol, mms, L): - xy = sol.nodes[:, :2] - return h1_semi_error( - xy, sol.conn, np.column_stack([sol.ux, sol.uy]), - lambda x, y: mms.grad_u_ex(x, y, L), - cls.ELEMENT_RULE) - - -class _QuadElement(_ElementBase2D): - LABEL = "Q1 quad" - ELEMENT_RULE = staticmethod(quad_q1_rule(2)) - - @staticmethod - def _source_rule(mms): - rule = mms.source_quadrature_quad - if rule is None: - raise ValueError( - f"{type(mms).__name__}.source_quadrature_quad must be set") - return rule - - @staticmethod - def add_topology(Beam): - topology = Beam.addObject("QuadSetTopologyContainer", name="topology", - quads="@../Grid/grid.quads", - position="@../Grid/grid.position") - Beam.addObject("QuadSetTopologyModifier") - return topology - - @staticmethod - def read_connectivity(topology): - return topology.quads.array().copy() - - @staticmethod - def to_triangles(conn): - tris = [] - for q in conn: - tris.append([q[0], q[1], q[2]]) - tris.append([q[0], q[2], q[3]]) - return np.array(tris) - - -class _TriElement(_ElementBase2D): - LABEL = "P1 tri" - ELEMENT_RULE = staticmethod(tri_p1_rule(3)) - - @staticmethod - def _source_rule(mms): - rule = mms.source_quadrature_tri - if rule is None: - raise ValueError( - f"{type(mms).__name__}.source_quadrature_tri must be set") - return rule - - @staticmethod - def add_topology(Beam): - topology = Beam.addObject("TriangleSetTopologyContainer", name="topology") - Beam.addObject("Quad2TriangleTopologicalMapping", - input="@../Grid/grid", output="@topology") - Beam.addObject("TriangleSetTopologyModifier") - return topology - - @staticmethod - def read_connectivity(topology): - return topology.triangles.array().copy() - - @staticmethod - def to_triangles(conn): - return conn - - -# --------------------------------------------------------------------------- -# 3D elements (Vec3d / full Hooke) -# --------------------------------------------------------------------------- - -class _ElementBase3D: - @classmethod - def compute_nodal_forces(cls, nodes_3d, conn, mms, L, E, nu, nx, ny, nz): - xyz = nodes_3d[:, :3] - - F = assemble_nodal_forces( - 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: - 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) - return F - - @classmethod - def compute_l2(cls, sol, mms, L): - return l2_error( - sol.nodes, sol.conn, - np.column_stack([sol.ux, sol.uy, sol.uz]), - lambda x, y, z: mms.u_ex(x, y, z, L), - cls.ELEMENT_RULE) - - @classmethod - def compute_h1(cls, sol, mms, L): - return h1_semi_error( - sol.nodes, sol.conn, - np.column_stack([sol.ux, sol.uy, sol.uz]), - lambda x, y, z: mms.grad_u_ex(x, y, z, L), - cls.ELEMENT_RULE) - - -class _HexElement(_ElementBase3D): - LABEL = "Q1 hex" - ELEMENT_RULE = staticmethod(hex_q1_rule(2)) - - @staticmethod - def _source_rule(mms): - rule = mms.source_quadrature_hex - if rule is None: - raise ValueError( - f"{type(mms).__name__}.source_quadrature_hex must be set") - return rule - - @staticmethod - def add_topology(Solid): - topology = Solid.addObject("HexahedronSetTopologyContainer", - name="topology", - hexahedra="@../Grid/grid.hexahedra", - position="@../Grid/grid.position") - Solid.addObject("HexahedronSetTopologyModifier") - return topology - - @staticmethod - def read_connectivity(topology): - return topology.hexahedra.array().copy() - - -# --------------------------------------------------------------------------- -# P1 tet — NOTE : compute_nodal_forces est SURCHARGÉE ici (pas héritée de -# _ElementBase3D) car la traction Neumann doit être assemblée sur des -# facettes TRIANGULAIRES (tri_face_rule / _boundary_triangles), cohérentes -# avec la base P1 réelle de l'élément — pas sur des facettes quad Q1 -# (quad_face_rule / _boundary_quads) comme le fait la version héritée. -# --------------------------------------------------------------------------- - -class _TetElement(_ElementBase3D): - LABEL = "P1 tet" - ELEMENT_RULE = staticmethod(tet_p1_rule(4)) - - @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 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() - - @classmethod - def compute_nodal_forces(cls, nodes_3d, conn, mms, L, E, nu, nx, ny, nz, - diagonal="main"): - xyz = nodes_3d[:, :3] - - F = assemble_nodal_forces( - lambda x, y, z: mms.source(x, y, z, E, nu, L), - xyz, conn, cls._source_rule(mms)) - - sides = _boundary_triangles_from_conn(conn, xyz, L) - xm, xp, ym, yp, zm, zp = (sides["xm"], sides["xp"], - sides["ym"], sides["yp"], - sides["zm"], sides["zp"]) - - tri_rule = tri_face_rule(3) - for tris, nrm_x, nrm_y, nrm_z in sides: - 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, tris, tri_rule) - return F - - - -def _boundary_triangles_from_conn(conn, xyz, L, eps=1e-9): - """Extrait les vraies faces de bord du maillage tet à partir de sa - connectivité (une face partagée par un seul tet = face de bord), - avec orientation sortante calculée géométriquement, puis les groupe - par face du cube. Contrairement à `_boundary_triangles`, ceci est - garanti cohérent avec la triangulation réelle produite par - Hexa2TetraTopologicalMapping, quel que soit `swapping`. - """ - face_local = [(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)] - count = {} - for tet in conn: - for loc in face_local: - key = tuple(sorted(tet[i] for i in loc)) - count[key] = count.get(key, 0) + 1 - - sides = {"xm": [], "xp": [], "ym": [], "yp": [], "zm": [], "zp": []} - for tet in conn: - for loc in face_local: - tri = [tet[i] for i in loc] - key = tuple(sorted(tri)) - if count[key] != 1: - continue - opp = [v for v in tet if v not in tri][0] - p0, p1, p2 = xyz[tri[0]], xyz[tri[1]], xyz[tri[2]] - n = np.cross(p1 - p0, p2 - p0) - if np.dot(n, xyz[opp] - p0) > 0: - tri = [tri[0], tri[2], tri[1]] - n = -n - nn = n / np.linalg.norm(n) - c = (xyz[tri[0]] + xyz[tri[1]] + xyz[tri[2]]) / 3.0 - if abs(c[0]) < eps and nn[0] < 0: sides["xm"].append(tuple(tri)) - elif abs(c[0] - L) < eps and nn[0] > 0: sides["xp"].append(tuple(tri)) - elif abs(c[1]) < eps and nn[1] < 0: sides["ym"].append(tuple(tri)) - elif abs(c[1] - L) < eps and nn[1] > 0: sides["yp"].append(tuple(tri)) - elif abs(c[2]) < eps and nn[2] < 0: sides["zm"].append(tuple(tri)) - elif abs(c[2] - L) < eps and nn[2] > 0: sides["zp"].append(tuple(tri)) - return sides - - -element_quad = _QuadElement() -element_tri = _TriElement() -element_hex = _HexElement() -element_tet = _TetElement() \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/fem.py b/examples/Freefem/MMS/3D_tetra/fem.py deleted file mode 100644 index 762a4bcc..00000000 --- a/examples/Freefem/MMS/3D_tetra/fem.py +++ /dev/null @@ -1,492 +0,0 @@ -"""FEM toolbox: quadrature, assembly, error norms. - -Two parallel families of routines: - -* **1D scalar-integrand** path used by the 1D bar driver: - `line_quadrature`, `assemble_nodal_forces_1d`, `l2_error_1d`, - `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), - `assemble_nodal_forces`, `assemble_traction`, - `l2_error`, `h1_semi_error`. - -Element-rule protocol: - - rule(xe) yields (coords, w, N, dN_phys) - xe : (n_local, dim) physical node coordinates of one element - coords : (dim,) physical Gauss-point coordinates - w : scalar weight (already × detJ / triangle area) - N : (n_local,) shape values at the Gauss point - dN_phys : (dim, n_local) physical-coord shape gradients, indexed - by (spatial axis, local node) - -Facet-rule protocol: - - rule(xe) yields (coords, w, N) - same convention; no shape gradient (not needed on boundary). -""" -import numpy as np - - -_GAUSS_LEGENDRE_1D = { - 1: (np.array([0.0]), - np.array([2.0])), - 2: (np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)]), - np.array([1.0, 1.0])), -} - - -# --------------------------------------------------------------------------- -# 1D scalar-integrand quadrature (consumed by the 1D bar driver only) -# --------------------------------------------------------------------------- - -def line_quadrature(n_pts): - """Return a Gauss-Legendre rule on a physical segment using `n_pts` points. - - The returned callable approximates int_{x1}^{x2} g(x) dx. - """ - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"line_quadrature: {n_pts}-point rule not supported") - xi, w = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(g, x1, x2): - h = x2 - x1 - x_k = 0.5 * (x1 + x2) + 0.5 * h * xi - return 0.5 * h * sum(w_i * g(x_i) for w_i, x_i in zip(w, x_k)) - return rule - - -# Error-norm rules. H1 requires >=2 pts: the P1 gradient superconverges at -# the element midpoint, so 1-pt fakes O(h^2) instead of O(h). -L2_QUADRATURE_1D = line_quadrature(2) -H1_QUADRATURE_1D = line_quadrature(2) - - -# --------------------------------------------------------------------------- -# Reference-element shape functions -# --------------------------------------------------------------------------- - -def _shape_q1(xi, eta): - """2D Q1 shape on [-1, 1]^2. Node order: (-,-), (+,-), (+,+), (-,+).""" - N = 0.25 * np.array([ - (1 - xi) * (1 - eta), - (1 + xi) * (1 - eta), - (1 + xi) * (1 + eta), - (1 - xi) * (1 + eta), - ]) - dN_dxi = 0.25 * np.array([-(1 - eta), (1 - eta), (1 + eta), -(1 + eta)]) - dN_deta = 0.25 * np.array([-(1 - xi), -(1 + xi), (1 + xi), (1 - xi) ]) - return N, dN_dxi, dN_deta - - -def _shape_hex_q1(xi, eta, zeta): - """3D Q1 hex shape on [-1, 1]^3. Node order matches SOFA RegularGridTopology: - 0:(-,-,-) 1:(+,-,-) 2:(+,+,-) 3:(-,+,-) - 4:(-,-,+) 5:(+,-,+) 6:(+,+,+) 7:(-,+,+).""" - N = 0.125 * np.array([ - (1 - xi) * (1 - eta) * (1 - zeta), - (1 + xi) * (1 - eta) * (1 - zeta), - (1 + xi) * (1 + eta) * (1 - zeta), - (1 - xi) * (1 + eta) * (1 - zeta), - (1 - xi) * (1 - eta) * (1 + zeta), - (1 + xi) * (1 - eta) * (1 + zeta), - (1 + xi) * (1 + eta) * (1 + zeta), - (1 - xi) * (1 + eta) * (1 + zeta), - ]) - dN_dxi = 0.125 * np.array([ - -(1 - eta) * (1 - zeta), (1 - eta) * (1 - zeta), - (1 + eta) * (1 - zeta), -(1 + eta) * (1 - zeta), - -(1 - eta) * (1 + zeta), (1 - eta) * (1 + zeta), - (1 + eta) * (1 + zeta), -(1 + eta) * (1 + zeta), - ]) - dN_deta = 0.125 * np.array([ - -(1 - xi) * (1 - zeta), -(1 + xi) * (1 - zeta), - (1 + xi) * (1 - zeta), (1 - xi) * (1 - zeta), - -(1 - xi) * (1 + zeta), -(1 + xi) * (1 + zeta), - (1 + xi) * (1 + zeta), (1 - xi) * (1 + zeta), - ]) - dN_dzeta = 0.125 * np.array([ - -(1 - xi) * (1 - eta), -(1 + xi) * (1 - eta), - -(1 + xi) * (1 + eta), -(1 - xi) * (1 + eta), - (1 - xi) * (1 - eta), (1 + xi) * (1 - eta), - (1 + xi) * (1 + eta), (1 - xi) * (1 + eta), - ]) - return N, dN_dxi, dN_deta, dN_dzeta - - -# Reference-triangle quadrature, integrating over {(xi,eta): xi,eta>=0, xi+eta<=1}. -# Weights sum to the reference area 1/2; the physical weight is w_ref * 2*area. -_TRI_QUADRATURE = { - 1: (np.array([[1/3, 1/3]]), - np.array([1/2])), - 3: (np.array([[1/6, 1/6], [2/3, 1/6], [1/6, 2/3]]), - np.array([1/6, 1/6, 1/6])), -} - -_a4, _b4 = 0.5854101966249685, 0.1381966011250105 # symmetric 4-point rule -_TET_QUADRATURE = { - 1: (np.array([[1/4, 1/4, 1/4]]), - np.array([1/6])), - 4: (np.array([[_a4, _b4, _b4], [_b4, _a4, _b4], - [_b4, _b4, _a4], [_b4, _b4, _b4]]), - np.array([1/24, 1/24, 1/24, 1/24])), -} - - -def quad_q1_rule(n_pts=2): - """Element rule for Q1 quads: tensor-product Gauss-Legendre on [-1,1]^2.""" - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"quad_q1_rule: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(xe): - xe_arr = np.asarray(xe, float) - x_col, y_col = xe_arr[:, 0], xe_arr[:, 1] - for xi, wi in zip(xi_pts, w_pts): - for eta, wj in zip(xi_pts, w_pts): - N, dN_dxi, dN_deta = _shape_q1(xi, eta) - J = np.array([[dN_dxi @ x_col, dN_dxi @ y_col], - [dN_deta @ x_col, dN_deta @ y_col]]) - detJ = np.linalg.det(J) - Jinv = np.linalg.inv(J) - coords = np.array([N @ x_col, N @ y_col]) - dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta - dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta - dN_phys = np.stack([dN_dx, dN_dy]) # (2, 4) [d, a] - yield coords, wi * wj * detJ, N, dN_phys - return rule - - -def tri_p1_rule(n_pts=1): - """Element rule for P1 triangles: 1-point (centroid) or 3-point Gauss. - - Shape gradients are constant per triangle. Node order in `xe` is the - canonical P1 ordering — corresponds to (N0, N1, N2) = (1-xi-eta, xi, eta). - """ - if n_pts not in _TRI_QUADRATURE: - raise ValueError(f"tri_p1_rule: {n_pts}-point rule not supported") - pts, wts = _TRI_QUADRATURE[n_pts] - - def rule(xe): - xe_arr = np.asarray(xe, float) # (3, 2) - (x0, y0), (x1, y1), (x2, y2) = xe_arr - A2 = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0) - area = abs(A2) / 2.0 - dN_dx = np.array([(y1 - y2) / A2, (y2 - y0) / A2, (y0 - y1) / A2]) - dN_dy = np.array([(x2 - x1) / A2, (x0 - x2) / A2, (x1 - x0) / A2]) - dN_phys = np.stack([dN_dx, dN_dy]) - for (xi, eta), w_ref in zip(pts, wts): - N = np.array([1.0 - xi - eta, xi, eta]) - coords = np.array([N @ xe_arr[:, 0], N @ xe_arr[:, 1]]) - yield coords, w_ref * 2.0 * area, N, dN_phys - return rule - - -def hex_q1_rule(n_pts=2): - """Element rule for Q1 hexes: tensor-product Gauss-Legendre on [-1,1]^3.""" - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"hex_q1_rule: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(xe): - xe_arr = np.asarray(xe, float) - x_col, y_col, z_col = xe_arr[:, 0], xe_arr[:, 1], xe_arr[:, 2] - for xi, wi in zip(xi_pts, w_pts): - for eta, wj in zip(xi_pts, w_pts): - for zeta, wk in zip(xi_pts, w_pts): - N, dN_dxi, dN_deta, dN_dzeta = _shape_hex_q1(xi, eta, zeta) - J = np.array([ - [dN_dxi @ x_col, dN_dxi @ y_col, dN_dxi @ z_col], - [dN_deta @ x_col, dN_deta @ y_col, dN_deta @ z_col], - [dN_dzeta @ x_col, dN_dzeta @ y_col, dN_dzeta @ z_col], - ]) - detJ = np.linalg.det(J) - Jinv = np.linalg.inv(J) - coords = np.array([N @ x_col, N @ y_col, N @ z_col]) - dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta + Jinv[2, 0] * dN_dzeta - dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta + Jinv[2, 1] * dN_dzeta - dN_dz = Jinv[0, 2] * dN_dxi + Jinv[1, 2] * dN_deta + Jinv[2, 2] * dN_dzeta - dN_phys = np.stack([dN_dx, dN_dy, dN_dz]) # (3, 8) [d, a] - yield coords, wi * wj * wk * detJ, N, dN_phys - return rule - - -def tet_p1_rule(n_pts=1): - """Element rule for P1 tetrahedra: 1-point (centroid) or 4-point Gauss. - - Shape gradients are constant per tetrahedron (linear element). Node - order in `xe` is the canonical P1 ordering — corresponds to - (N0, N1, N2, N3) = (1-xi-eta-zeta, xi, eta, zeta), matching SOFA's - `Hexa2TetraTopologicalMapping` tetrahedron orientation. - """ - 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] - - def rule(xe): - xe_arr = np.asarray(xe, float) # (4, 3) - x0, x1, x2, x3 = xe_arr - # Constant Jacobian: edge vectors from node 0 as columns. - J = np.column_stack([x1 - x0, x2 - x0, x3 - x0]) # (3, 3) - detJ = np.linalg.det(J) - vol = abs(detJ) / 6.0 - Jinv = np.linalg.inv(J) - # Reference gradients of (N0, N1, N2, N3) w.r.t. (xi, eta, zeta). - dN_dxi = np.array([-1.0, 1.0, 0.0, 0.0]) - dN_deta = np.array([-1.0, 0.0, 1.0, 0.0]) - dN_dzeta = np.array([-1.0, 0.0, 0.0, 1.0]) - dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta + Jinv[2, 0] * dN_dzeta - dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta + Jinv[2, 1] * dN_dzeta - dN_dz = Jinv[0, 2] * dN_dxi + Jinv[1, 2] * dN_deta + Jinv[2, 2] * dN_dzeta - dN_phys = np.stack([dN_dx, dN_dy, dN_dz]) # (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) -# --------------------------------------------------------------------------- - -def edge_line_rule(n_pts=2): - """Facet rule for a 2D boundary edge: linear shape, 1D Gauss-Legendre. - - Yields per Gauss point along the edge: (coords (2,), w, N (2,)). - """ - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"edge_line_rule: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(xe): - xe_arr = np.asarray(xe, float) # (2, 2) - Le = np.linalg.norm(xe_arr[1] - xe_arr[0]) - for xi, wi in zip(xi_pts, w_pts): - t = 0.5 * (xi + 1.0) - N = np.array([1.0 - t, t]) - coords = N @ xe_arr # (2,) - yield coords, wi * Le / 2.0, N - return rule - - -def quad_face_rule(n_pts=4): - """Facet rule for a 3D boundary quad face: bilinear shape, 2D Gauss. - - Yields per Gauss point on the face: (coords (3,), w, N (4,)). - Quad orientation does not affect the magnitude `dA = ||t_xi × t_eta||`. - """ - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"quad_face_rule: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(xe): - xe_arr = np.asarray(xe, float) # (4, 3) - for xi, wi in zip(xi_pts, w_pts): - for eta, wj in zip(xi_pts, w_pts): - N = 0.25 * np.array([ - (1 - xi) * (1 - eta), - (1 + xi) * (1 - eta), - (1 + xi) * (1 + eta), - (1 - xi) * (1 + eta), - ]) - dN_dxi = 0.25 * np.array([-(1 - eta), (1 - eta), (1 + eta), -(1 + eta)]) - dN_deta = 0.25 * np.array([-(1 - xi), -(1 + xi), (1 + xi), (1 - xi) ]) - t_xi = dN_dxi @ xe_arr # (3,) - t_eta = dN_deta @ xe_arr - dA = np.linalg.norm(np.cross(t_xi, t_eta)) - coords = N @ xe_arr # (3,) - yield coords, wi * wj * dA, N - return rule - - - -def tri_face_rule(n_pts=3): - """Facet rule for a 3D boundary triangle: linear shape, 2D triangle Gauss. - - Yields per Gauss point on the face: (coords (3,), w, N (3,)). - """ - 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) - v1 = xe_arr[1] - xe_arr[0] - v2 = xe_arr[2] - xe_arr[0] - area2 = np.linalg.norm(np.cross(v1, v2)) # = 2 * aire du triangle - for (xi, eta), w_ref in zip(pts, wts): - N = np.array([1.0 - xi - eta, xi, eta]) - coords = N @ xe_arr - yield coords, w_ref * area2, N - return rule - - -# --------------------------------------------------------------------------- -# 1D FEM assembly (scalar-integrand path) -# --------------------------------------------------------------------------- - -def assemble_nodal_forces_1d(f_body, nodes, edges, quadrature): - """ - Assemble the consistent nodal force vector F_i = integral f_body(x) phi_i(x) dx. - - f_body : callable x -> float - nodes : 1-D array of node coordinates - edges : iterable of (a, b) node-index pairs defining each element - quadrature: rule from line_quadrature(n_pts) - """ - forces = np.zeros(len(nodes)) - for a, b in edges: - x1, x2 = nodes[a], nodes[b] - h = x2 - x1 - forces[a] += quadrature(lambda x, x1=x1, x2=x2, h=h: f_body(x) * (x2 - x) / h, x1, x2) - forces[b] += quadrature(lambda x, x1=x1, x2=x2, h=h: f_body(x) * (x - x1) / h, x1, x2) - return forces - - -# --------------------------------------------------------------------------- -# Dim-agnostic FEM assembly (2D and 3D) -# --------------------------------------------------------------------------- - -def assemble_nodal_forces(f_body, nodes, conn, element_rule): - """ - Assemble consistent nodal body forces, dim- and element-type agnostic. - - F_a = sum_e integral_{Omega_e} N_a(x) * f_body(x) dx - - f_body : callable (*coords) -> dim-tuple of force components - nodes : (N, dim) array - conn : iterable of node-index lists, one per element - element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule - """ - nodes = np.asarray(nodes) - dim = nodes.shape[1] - F = np.zeros((len(nodes), dim)) - for elem in conn: - - idx = np.asarray(elem) - xe = nodes[idx] # (n_local, dim) - for coords, w, N, _ in element_rule(xe): - f_components = f_body(*coords) - for a, node in enumerate(idx): - for d in range(dim): - F[node, d] += N[a] * f_components[d] * w - return F - - -def assemble_traction(traction, nodes, facets, facet_rule): - """ - Assemble consistent nodal forces from a boundary traction over facets. - - F_a += integral_{facet} N_a(s) * traction(x(s)) dS - - The outward unit normal is baked into `traction` by the caller (one - lambda per boundary face); the facet rule only supplies the surface area - element and shape values. - - traction : callable (*coords) -> dim-tuple of force components - nodes : (N, dim) array - facets : iterable of node-index tuples (one per boundary facet) - facet_rule : rule from edge_line_rule (2D) / quad_face_rule (3D) - """ - nodes = np.asarray(nodes) - dim = nodes.shape[1] - F = np.zeros((len(nodes), dim)) - for facet in facets: - idx = np.asarray(facet) - xe = nodes[idx] # (n_local, dim) - for coords, w, N in facet_rule(xe): - T = traction(*coords) - for a, node in enumerate(idx): - for d in range(dim): - F[node, d] += N[a] * T[d] * w - return F - - -# --------------------------------------------------------------------------- -# 1D error norms -# --------------------------------------------------------------------------- - -def l2_error_1d(nodes, edges, u_h, u_ex, quadrature): - """L2 error norm: sqrt( integral (u_h - u_ex)^2 dx ) over the mesh.""" - total = 0.0 - for a, b in edges: - x1, x2 = nodes[a], nodes[b] - h = x2 - x1 - u_a, u_b = u_h[a], u_h[b] - u_interp = lambda x, x1=x1, h=h, u_a=u_a, u_b=u_b: u_a + (u_b - u_a) * (x - x1) / h - total += quadrature(lambda x: (u_interp(x) - u_ex(x)) ** 2, x1, x2) - return np.sqrt(total) - - -def h1_semi_error_1d(nodes, edges, u_h, du_ex, quadrature): - """H1 semi-norm error: sqrt( integral (du_h - du_ex)^2 dx ) over the mesh.""" - total = 0.0 - for a, b in edges: - x1, x2 = nodes[a], nodes[b] - h = x2 - x1 - # du_h = sum_a u_a * dphi_a/dx, with dphi_a/dx = -1/h, dphi_b/dx = +1/h - du_h = u_h[a] * (-1.0 / h) + u_h[b] * (1.0 / h) - total += quadrature(lambda x, du_h=du_h: (du_h - du_ex(x)) ** 2, x1, x2) - return np.sqrt(total) - - -# --------------------------------------------------------------------------- -# Dim-agnostic error norms (2D and 3D) -# --------------------------------------------------------------------------- - -def l2_error(nodes, conn, u_h, u_ex, element_rule): - """ - Vector L2 error norm, dim- and element-type agnostic. - - sqrt( integral || u_h(x) - u_ex(x) ||^2 dx ) over the mesh - - nodes : (N, dim) array - conn : iterable of node-index lists, one per element - u_h : (N, dim) array of nodal displacements - u_ex : callable (*coords) -> dim-tuple - element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule - """ - nodes = np.asarray(nodes) - u = np.asarray(u_h) - total = 0.0 - for elem in conn: - idx = np.asarray(elem) - xe = nodes[idx] # (n_local, dim) - u_loc = u[idx] # (n_local, dim) - for coords, w, N, _ in element_rule(xe): - u_h_g = N @ u_loc # (dim,) - u_ex_g = np.asarray(u_ex(*coords)) # (dim,) - diff = u_h_g - u_ex_g - total += float(np.dot(diff, diff)) * w - return np.sqrt(total) - - -def h1_semi_error(nodes, conn, u_h, grad_u_ex, element_rule): - """ - Vector H1 semi-norm error, dim- and element-type agnostic. - - sqrt( integral || grad u_h - grad u_ex ||_F^2 dx ) over the mesh - - nodes : (N, dim) array - conn : iterable of node-index lists, one per element - u_h : (N, dim) array of nodal displacements - grad_u_ex : callable (*coords) -> (dim, dim) array with [i, j] = du_i/dx_j - element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule - """ - nodes = np.asarray(nodes) - u = np.asarray(u_h) - total = 0.0 - for elem in conn: - idx = np.asarray(elem) - xe = nodes[idx] # (n_local, dim) - u_loc = u[idx] # (n_local, dim) - for coords, w, _, dN_phys in element_rule(xe): - # grad_uh[i, d] = sum_a u_i[a] * dN_a/dx_d - grad_uh = (dN_phys @ u_loc).T # (dim, dim) [i, d] - grad_ue = np.asarray(grad_u_ex(*coords)) # [i, j] = du_i/dx_j - diff = grad_uh - grad_ue - total += float(np.sum(diff * diff)) * w - return np.sqrt(total) \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/manufactured_solution.py b/examples/Freefem/MMS/3D_tetra/manufactured_solution.py deleted file mode 100644 index 7e6c603d..00000000 --- a/examples/Freefem/MMS/3D_tetra/manufactured_solution.py +++ /dev/null @@ -1,63 +0,0 @@ -"""Abstract base class for 3D MMS cases (Cartesian domain [0,L]^3).""" - -import os -import sys -from abc import abstractmethod -import numpy as np - -sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) -from mms_case import MMSCase - - -def lame(E, nu): - """Lamé parameters (lambda, mu) for 3D linear elasticity (full Hooke).""" - lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) - mu = E / (2.0 * (1.0 + nu)) - return lam, mu - - -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): - """Exact solution: returns (ux, uy, uz).""" - - @abstractmethod - def grad_u_ex(self, x, y, z, L): - """Exact 3x3 displacement gradient - [[dux/dx, dux/dy, dux/dz], - [duy/dx, duy/dy, duy/dz], - [duz/dx, duz/dy, duz/dz]].""" - - @abstractmethod - def source(self, x, y, z, E, nu, L): - """Body force: returns (fx, fy, fz).""" - - @abstractmethod - def apply_bcs(self, Solid, nodes_3d, L): - """Install Dirichlet constraints on the SOFA `Solid` node. - - The MMS author chooses the BC pattern matching the manufactured - solution. Neumann tractions on the six faces are already assembled - into the consistent nodal force by the framework (via `self.traction`); - this method only needs to add SOFA constraint objects. - """ - - def traction(self, x, y, z, nx, ny, nz, E, nu, L): - """Traction on a face with outward unit normal (nx, ny, nz): - returns (Tx, Ty, Tz). - - Derived from `grad_u_ex` via sigma = lambda tr(eps) I + 2 mu eps and - T = sigma . n. - """ - lam, mu = lame(E, nu) - G = self.grad_u_ex(x, y, z, L) - eps = 0.5 * (G + G.T) - tr = eps[0, 0] + eps[1, 1] + eps[2, 2] - S = lam * tr * np.eye(3) + 2.0 * mu * eps - T = S @ np.array([nx, ny, nz]) - return T[0], T[1], T[2] \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/mms_case.py b/examples/Freefem/MMS/3D_tetra/mms_case.py deleted file mode 100644 index 8f426ace..00000000 --- a/examples/Freefem/MMS/3D_tetra/mms_case.py +++ /dev/null @@ -1,15 +0,0 @@ -"""Cross-dim base class for manufactured-solution cases. - -Dim-specific bases (`MMSCase1D`, `MMSCase2D`, `MMSCase3D`) inherit from -`MMSCase` and declare the abstract solution/derivative/source/BC methods -plus their `source_quadrature_*` placeholders. The case files -(`cubic.py`, `trigonometric.py`, `sinus_neumann.py`, …) inherit from the -dim-specific bases and override the abstract methods explicitly. -""" - -from abc import ABC - - -class MMSCase(ABC): - name = None # case identifier (must match the params.json key) - plot_label = None # LaTeX label for the exact solution diff --git a/examples/Freefem/MMS/3D_tetra/output.py b/examples/Freefem/MMS/3D_tetra/output.py deleted file mode 100644 index 29d7c34c..00000000 --- a/examples/Freefem/MMS/3D_tetra/output.py +++ /dev/null @@ -1,105 +0,0 @@ -"""Shared output helpers across MMS 1D/2D/3D drivers (tables + plots).""" - -import os -import numpy as np -import matplotlib.pyplot as plt - -from plot_utils import annotate_convergence_rates - - -_AXES = ("x", "y", "z") - - -def write_solution_table(stem, coords, u_h, u_ex, results_dir, error_dict): - """ - Write per-node solution table and error summary to /.txt. - - coords : (N,) or (N, dim) array of node coordinates - u_h : (N,) or (N, dim) array of nodal displacements - u_ex : callable (*coords_row) -> scalar (dim=1) or dim-tuple - error_dict : ordered dict of {label: value} for error summary lines - - Column count adjusts with dim. For 1D the displacement column header is - `ux_h` and the error header is `err_x` (uniform with 2D/3D). - """ - coords = np.asarray(coords) - if coords.ndim == 1: - coords = coords.reshape(-1, 1) - u_h = np.asarray(u_h) - if u_h.ndim == 1: - u_h = u_h.reshape(-1, 1) - dim = coords.shape[1] - assert u_h.shape[1] == dim, "u_h column count must match coord dim" - axes = _AXES[:dim] - - coord_hdr = " | ".join(f"{ax:>10}" for ax in axes) - disp_hdr = " | ".join(f"{'u'+ax+'_h':>15}" for ax in axes) - exact_hdr = " | ".join(f"{'u'+ax+'_ex':>15}" for ax in axes) - err_hdr = " | ".join(f"{'err_'+ax:>15}" for ax in axes) - header = f"{coord_hdr} | {disp_hdr} | {exact_hdr} | {err_hdr}" - - os.makedirs(results_dir, exist_ok=True) - path = os.path.join(results_dir, f"{stem}.txt") - with open(path, "w") as f: - f.write(header + "\n") - f.write("-" * len(header) + "\n") - for c_row, uh_row in zip(coords, u_h): - ue = u_ex(*c_row) - ue = (ue,) if np.isscalar(ue) else tuple(ue) - coord_str = " | ".join(f"{x:10.4f}" for x in c_row) - disp_str = " | ".join(f"{x:15.6e}" for x in uh_row) - exact_str = " | ".join(f"{x:15.6e}" for x in ue) - err_str = " | ".join(f"{abs(uh_row[k] - ue[k]):15.6e}" - for k in range(dim)) - f.write(f"{coord_str} | {disp_str} | {exact_str} | {err_str}\n") - f.write("\n") - for label, val in error_dict.items(): - f.write(f"{label:12s} = {val:.6e}\n") - - -def write_convergence_table(stem, rows, results_dir): - """ - Write convergence table to /.txt. - - rows : list of dicts with keys 'nx', 'h', and one key per error column. - Rate columns are strings (empty for the first row). - """ - os.makedirs(results_dir, exist_ok=True) - path = os.path.join(results_dir, f"{stem}.txt") - err_keys = [k for k in rows[0] if k not in ("nx", "h")] - header = f"{'nx':>6} | {'h':>10}" + "".join(f" | {k:>16}" for k in err_keys) - with open(path, "w") as f: - f.write(header + "\n") - f.write("-" * len(header) + "\n") - for row in rows: - line = f"{row['nx']:6d} | {row['h']:10.4f}" - for k in err_keys: - v = row[k] - line += f" | {v:16.6e}" if isinstance(v, float) else f" | {v:>16}" - f.write(line + "\n") - - -def plot_convergence(stem, hs, series, title, results_dir, ylabel="Error"): - """ - Save log-log convergence plot to /.png. - - series : list of {"label", "errors", "style"?} dicts - Per-segment convergence rates are annotated above each line segment. - """ - os.makedirs(results_dir, exist_ok=True) - h_arr = np.array(hs) - default = ["bo-", "rs--", "g^:", "m^-"] - fig, ax = plt.subplots(figsize=(8, 5)) - for i, s in enumerate(series): - style = s.get("style", default[i % len(default)]) - e_arr = np.array(s["errors"]) - ax.loglog(h_arr, e_arr, style, label=s["label"], linewidth=2, markersize=7) - annotate_convergence_rates(ax, h_arr, e_arr) - ax.set_xlabel("h") - ax.set_ylabel(ylabel) - ax.set_title(title) - ax.legend() - ax.grid(True, alpha=0.3, which="both") - fig.tight_layout() - fig.savefig(os.path.join(results_dir, f"{stem}.png"), dpi=150) - plt.close(fig) diff --git a/examples/Freefem/MMS/3D_tetra/params.json b/examples/Freefem/MMS/3D_tetra/params.json deleted file mode 100644 index 7d0b23b4..00000000 --- a/examples/Freefem/MMS/3D_tetra/params.json +++ /dev/null @@ -1,19 +0,0 @@ -{ - "length": 1.0, - "youngModulus": 1e6, - "forceField": "LinearSmallStrainFEMForceField", - "linearSolver": { - "type": "PCGLinearSolver", - "parameters": { "iterations": 5000, "tolerance": 1e-14 } - }, - "reference": { - "nx": 20, - "nu": 0.3 - }, - "convergence": { - "nu_values": [0.3], - "nx_values": { - "sinusoidal": [5,10,12,16,20,25,50] - } - } -} \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/plot_utils.py b/examples/Freefem/MMS/3D_tetra/plot_utils.py deleted file mode 100644 index 5764c32d..00000000 --- a/examples/Freefem/MMS/3D_tetra/plot_utils.py +++ /dev/null @@ -1,16 +0,0 @@ -"""Plot helpers shared across MMS 1D/2D drivers.""" - -import numpy as np - - -def annotate_convergence_rates(ax, hs, errors): - """Annotate the per-segment log-log convergence rate above each segment.""" - hs = np.asarray(hs) - errors = np.asarray(errors) - for k in range(1, len(hs)): - rate = np.log(errors[k] / errors[k-1]) / np.log(hs[k] / hs[k-1]) - x_mid = np.sqrt(hs[k] * hs[k-1]) - y_mid = np.sqrt(errors[k] * errors[k-1]) - ax.annotate(f"{rate:.2f}", xy=(x_mid, y_mid), - xytext=(0, 8), textcoords="offset points", - ha="center", fontsize=9) diff --git a/examples/Freefem/MMS/3D_tetra/run_convergence.py b/examples/Freefem/MMS/3D_tetra/run_convergence.py deleted file mode 100644 index ea27e5ef..00000000 --- a/examples/Freefem/MMS/3D_tetra/run_convergence.py +++ /dev/null @@ -1,85 +0,0 @@ -"""Run the mesh-refinement convergence study for every 3D MMS case. - -Loops cases × nu × nx and writes per-(case, nu) text tables and convergence -plots into the shared `results/` directory. Mirrors the 2D driver minus the -plane-stress / plane-strain `dim` axis (3D has a single constitutive branch). -""" - -from sinusoidal import mms as sinusoidal_mms - -from solid import ( - RESULTS_DIR, - load_params, - solve_solid, - element_hex, - element_tet, -) -from convergence import run_convergence_series -from output import plot_convergence - - -def convergence_study(elem_specs, mms, L, E, nu, nx_values, - force_field, linear_solver): - """ - Run a convergence series for each element type in elem_specs, write a - per-(element) text table, and one shared plot with L²/H¹ for every - element on the same axes. - - elem_specs : list of dicts with keys 'elem', 'label', 'l2_style', 'h1_style' - force_field : name of the FEM force field to test - """ - print(f"\n PoissonRatio = {nu}", flush=True) - - plot_series, hs_ref = [], None - for spec in elem_specs: - elem, label = spec["elem"], spec["label"] - tag = label.replace(" ", "_") - stem = f"convergence_{mms.name}_{tag}_nu{nu}" - - hs, errors = run_convergence_series( - nx_values = nx_values, - run_fn = lambda nx, _e=elem: solve_solid( - _e, mms, L, E, nu, nx, nx, nx, - force_field=force_field, linear_solver=linear_solver), - h_fn = lambda nx: L / (nx - 1), - error_fns = { - "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L), - "H1": lambda sol, _e=elem: _e.compute_h1(sol, mms, L), - }, - banner = f"-- {label} {mms.name} nu={nu} --", - results_dir = RESULTS_DIR, - table_stem = stem, - ) - - plot_series.append({"label": f"{label} L²", - "errors": errors["L2"], "style": spec["l2_style"]}) - plot_series.append({"label": f"{label} H¹", - "errors": errors["H1"], "style": spec["h1_style"]}) - hs_ref = hs - - title = f"Convergence — {mms.name} nu={nu}" - plot_convergence(f"convergence_{mms.name}_nu{nu}", - hs_ref, plot_series, title=title, results_dir=RESULTS_DIR) - - -if __name__ == "__main__": - cfg = load_params() - L = cfg["length"] - E = cfg["youngModulus"] - ff = cfg["forceField"] - ls = cfg["linearSolver"] - conv = cfg["convergence"] - - specs = [ - {"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,): - nx_vals = conv["nx_values"][mms.name] - print(f"\n== {mms.name} ==") - for nu in conv["nu_values"]: - convergence_study(specs, mms, L, E, nu, nx_vals, - force_field=ff, linear_solver=ls) \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/run_convergence_dirichlet.py b/examples/Freefem/MMS/3D_tetra/run_convergence_dirichlet.py deleted file mode 100644 index a15e36e1..00000000 --- a/examples/Freefem/MMS/3D_tetra/run_convergence_dirichlet.py +++ /dev/null @@ -1,40 +0,0 @@ -"""Diagnostic : convergence P1 tet avec Dirichlet complet (isole le Neumann).""" - -from convergence import run_convergence_series -from output import plot_convergence -from solid import load_params, RESULTS_DIR, element_tet -from sinusoidal import mms as sinusoidal_mms -from sinusoidal_dirichlet import solve_dirichlet - - -if __name__ == "__main__": - cfg = load_params() - L, E = cfg["length"], cfg["youngModulus"] - ff = cfg["forceField"] - ls = cfg["linearSolver"] - conv = cfg["convergence"] - - nx_vals = conv["nx_values"][sinusoidal_mms.name] - nu = conv["nu_values"][0] - - hs, errors = run_convergence_series( - nx_values = nx_vals, - run_fn = lambda nx: solve_dirichlet( - element_tet, sinusoidal_mms, L, E, nu, nx, nx, nx, - force_field=ff, linear_solver=ls), - h_fn = lambda nx: L / (nx - 1), - error_fns = { - "L2": lambda sol: element_tet.compute_l2(sol, sinusoidal_mms, L), - "H1": lambda sol: element_tet.compute_h1(sol, sinusoidal_mms, L), - }, - banner = f"-- P1 tet DIRICHLET-only {sinusoidal_mms.name} nu={nu} --", - results_dir = RESULTS_DIR, - table_stem = f"convergence_{sinusoidal_mms.name}_P1tet_DIRICHLET_nu{nu}", - ) - - plot_convergence( - f"convergence_{sinusoidal_mms.name}_DIRICHLET_nu{nu}", - hs, [{"label": "P1 tet L² (Dirichlet)", "errors": errors["L2"], "style": "g^-"}, - {"label": "P1 tet H¹ (Dirichlet)", "errors": errors["H1"], "style": "md--"}], - title=f"Convergence — {sinusoidal_mms.name} DIRICHLET nu={nu}", - results_dir=RESULTS_DIR) \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/scene.py b/examples/Freefem/MMS/3D_tetra/scene.py deleted file mode 100644 index bcd390fd..00000000 --- a/examples/Freefem/MMS/3D_tetra/scene.py +++ /dev/null @@ -1,40 +0,0 @@ -"""Shared SOFA scene helpers (controllers) across MMS 1D/2D/3D drivers.""" - -import Sofa -import Sofa.Core - - -class NodalForceAssembler(Sofa.Core.Controller): - """Fill a placeholder ConstantForceField after Sofa init has run. - - SOFA topology components (RegularGridTopology, - EdgeSetTopologyContainer, QuadSetTopologyContainer, - HexahedronSetTopologyContainer, …) are only populated after init runs, - not during Python scene-build time. This controller defers nodal-force - assembly to onSimulationInitDoneEvent, where it reads rest positions - off the MechanicalObject and hands them to `compute_forces` together - with the topology container. The returned array is written into the - force field. - - Parameters - ---------- - dofs : SOFA MechanicalObject with rest_position - topology : SOFA topology container (read by compute_forces) - force_field : SOFA ConstantForceField to populate - compute_forces : callable (nodes, topology) -> ndarray sized like - force_field.forces - """ - - def __init__(self, dofs, topology, force_field, compute_forces, - *args, **kwargs): - super().__init__(*args, **kwargs) - self.dofs = dofs - self.topology = topology - self.force_field = force_field - self.compute_forces = compute_forces - - def onSimulationInitDoneEvent(self, event): - nodes = self.dofs.rest_position.array().copy() - F = self.compute_forces(nodes, self.topology) - with self.force_field.forces.writeableArray() as forces: - forces[:] = F diff --git a/examples/Freefem/MMS/3D_tetra/sinusoidal.py b/examples/Freefem/MMS/3D_tetra/sinusoidal.py deleted file mode 100644 index 00a1e9bb..00000000 --- a/examples/Freefem/MMS/3D_tetra/sinusoidal.py +++ /dev/null @@ -1,113 +0,0 @@ -""" -Sinusoidal 3D MMS on [0,L]^3 with linear-elasticity constitutive law: - - u_ex(x, y, z) = A * ( sin(pi x / L) sin(pi y / L), - sin(pi y / L) sin(pi z / L), - sin(pi z / L) sin(pi x / L) ) - - sigma = lambda tr(eps) I + 2 mu eps (3D / full Hooke). - -BCs: minimal Dirichlet that removes the six rigid-body modes — - - corner (0,0,0): pin all three components. - - corner (L,0,0): pin u_y and u_z (kill rotation about z and y). - - corner (0,L,0): pin u_z (kill rotation about x). -Everything else is Neumann on the six faces (assembled by the framework -from `self.traction`). -""" - -import numpy as np - -from manufactured_solution import MMSCase3D, lame -from solid import (case_scene, run_reference_scene, - element_hex, element_tet, hex_q1_rule) -from fem import tet_p1_rule - - -SINUS_AMPLITUDE = 1e-3 - - -class SinusNeumann(MMSCase3D): - name = "sinusoidal" - plot_label = (r"$u = A(\sin(\pi x/L)\sin(\pi y/L),\ " - r"\sin(\pi y/L)\sin(\pi z/L),\ " - 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 - A = SINUS_AMPLITUDE - return (A * np.sin(k*x) * np.sin(k*y), - A * np.sin(k*y) * np.sin(k*z), - A * np.sin(k*z) * np.sin(k*x)) - - def grad_u_ex(self, x, y, z, L): - k = np.pi / L - A = SINUS_AMPLITUDE - zero = 0.0 if np.isscalar(x) else np.zeros_like(np.asarray(x, float)) - dux_dx = A * k * np.cos(k*x) * np.sin(k*y) - dux_dy = A * k * np.sin(k*x) * np.cos(k*y) - dux_dz = zero - duy_dx = zero - duy_dy = A * k * np.cos(k*y) * np.sin(k*z) - duy_dz = A * k * np.sin(k*y) * np.cos(k*z) - duz_dx = A * k * np.cos(k*x) * np.sin(k*z) - duz_dy = zero - duz_dz = A * k * np.cos(k*z) * np.sin(k*x) - return np.array([[dux_dx, dux_dy, dux_dz], - [duy_dx, duy_dy, duy_dz], - [duz_dx, duz_dy, duz_dz]]) - - def source(self, x, y, z, E, nu, L): - lam, mu = lame(E, nu) - A = SINUS_AMPLITUDE - p = np.pi / L - sx, sy, sz = np.sin(p*x), np.sin(p*y), np.sin(p*z) - cx, cy, cz = np.cos(p*x), np.cos(p*y), np.cos(p*z) - - # Component-wise Laplacian - lap_ux = -2 * p**2 * sx * sy - lap_uy = -2 * p**2 * sy * sz - lap_uz = -2 * p**2 * sz * sx - - # Gradient of the divergence - d_divu_dx = p**2 * (-sx*sy + cz*cx) - d_divu_dy = p**2 * ( cx*cy - sy*sz) - d_divu_dz = p**2 * ( cy*cz - sz*sx) - - fx = A * (-(lam + mu) * d_divu_dx - mu * lap_ux) - fy = A * (-(lam + mu) * d_divu_dy - mu * lap_uy) - fz = A * (-(lam + mu) * d_divu_dz - mu * lap_uz) - return (fx, fy, fz) - - def apply_bcs(self, Solid, nodes_3d, L): - eps = 1e-10 - xyz = nodes_3d[:, :3] - - 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") - - 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) - - 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") - - -mms = SinusNeumann() -createScene = case_scene(mms,element_tet) - - -if __name__ == "__main__": - run_reference_scene(element_tet, mms) \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/sinusoidal_dirichlet.py b/examples/Freefem/MMS/3D_tetra/sinusoidal_dirichlet.py deleted file mode 100644 index 2487c04f..00000000 --- a/examples/Freefem/MMS/3D_tetra/sinusoidal_dirichlet.py +++ /dev/null @@ -1,240 +0,0 @@ -import numpy as np -import Sofa -import Sofa.Core -import Sofa.Simulation - -from solid import (get_nodes_3d, element_tet, load_params, - plot_solution_profile, plot_solution_slices, RESULTS_DIR) -from solid_solution import SolidSolution3D -from scene import NodalForceAssembler -from fem import assemble_nodal_forces -from output import write_solution_table -from sinusoidal import mms # réutilise l'instance SinusNeumann déjà définie là-bas - - -# --------------------------------------------------------------------------- -# Boundary helpers -# --------------------------------------------------------------------------- - -def _boundary_node_indices(nx, ny, nz): - """Tous les indices de nœuds sur une des 6 faces de la grille structurée.""" - def idx(i, j, k): - return i + j * nx + k * nx * ny - - idxs = set() - for j in range(ny): - for k in range(nz): - idxs.add(idx(0, j, k)); idxs.add(idx(nx - 1, j, k)) - for i in range(nx): - for k in range(nz): - idxs.add(idx(i, 0, k)); idxs.add(idx(i, ny - 1, k)) - for i in range(nx): - for j in range(ny): - idxs.add(idx(i, j, 0)); idxs.add(idx(i, j, nz - 1)) - return sorted(idxs) - - -def _body_force_only(element, mms, L, E, nu, nx, ny, nz): - """Forces nodales : uniquement le terme source volumique, pas de traction.""" - def compute(nodes, topology): - conn = element.read_connectivity(topology) - xyz = nodes[:, :3] - return assemble_nodal_forces( - lambda x, y, z: mms.source(x, y, z, E, nu, L), - xyz, conn, element._source_rule(mms)) - return compute - - -# --------------------------------------------------------------------------- -# Controller: pushes the Dirichlet boundary positions at first animate step -# --------------------------------------------------------------------------- - -class DirichletApplier(Sofa.Core.Controller): - """Ecrit u_ex sur les DOFs de bord au premier pas d'animation. - - La scène charge donc plate (position == rest_position partout au GUI - init), et la déformation de bord n'apparaît qu'au clic sur Play/Step, - juste avant que le NewtonRaphsonSolver ne résolve l'intérieur. - """ - - def __init__(self, dofs, boundary_idx, boundary_pos, *args, **kwargs): - super().__init__(*args, **kwargs) - self.dofs = dofs - self.boundary_idx = boundary_idx - self.boundary_pos = boundary_pos - self._applied = False - - def onAnimateBeginEvent(self, event): - if self._applied: - return - with self.dofs.position.writeable() as pos: - for i, p in zip(self.boundary_idx, self.boundary_pos): - pos[i] = p - self._applied = True - - -# --------------------------------------------------------------------------- -# Scene builder -# --------------------------------------------------------------------------- - -def build_dirichlet_scene(rootNode, mms, element, force_field, linear_solver, - L=1.0, E=1e6, nu=0.3, nx=6, ny=6, nz=6, - with_visual=False): - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", - "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", - "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.LinearSolver.Iterative", - "Sofa.Component.MechanicalLoad", - "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Grid", - "Sofa.Component.Topology.Container.Dynamic", - "Sofa.Component.Topology.Mapping", - "Sofa.Component.Visual", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_3d = get_nodes_3d(L, nx, ny, nz) - - Grid = rootNode.addChild("Grid") - Grid.addObject("RegularGridTopology", name="grid", - nx=nx, ny=ny, nz=nz, - min=[0.0, 0.0, 0.0], max=[L, L, L]) - - Solid = rootNode.addChild("Solid") - Solid.addObject("StaticSolver", name="staticSolver", printLog=False) - Solid.addObject("NewtonRaphsonSolver", name="newtonSolver", - maxNbIterationsNewton=1, - absoluteResidualStoppingThreshold=1e-10, - printLog=False) - Solid.addObject(linear_solver["type"], name="linearSolver", - **linear_solver["parameters"]) - - # Position/valeur cible du Dirichlet non-homogène sur le bord. - boundary_idx = _boundary_node_indices(nx, ny, nz) - init_pos = nodes_3d.copy() - for i in boundary_idx: - x, y, z = nodes_3d[i] - ux, uy, uz = mms.u_ex(x, y, z, L) - init_pos[i] = [x + ux, y + uy, z + uz] - boundary_pos = [init_pos[i].tolist() for i in boundary_idx] - - # MechanicalObject démarre à rest partout (rendu "plat" au chargement). - 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") - - Solid.addObject("FixedProjectiveConstraint", name="fix_boundary", - indices=boundary_idx) - - n_nodes = len(nodes_3d) - ff = Solid.addObject("ConstantForceField", name="MMS_forces", - template="Vec3d", - indices=list(range(n_nodes)), - forces=[[0.0, 0.0, 0.0]] * n_nodes) - - Solid.addObject(NodalForceAssembler( - dofs=dofs, topology=topology, force_field=ff, - compute_forces=_body_force_only(element, mms, L, E, nu, nx, ny, nz), - name="nodalForceAssembler")) - - Solid.addObject(DirichletApplier( - dofs=dofs, boundary_idx=boundary_idx, boundary_pos=boundary_pos, - name="dirichletApplier")) - - return dofs, topology - - -# --------------------------------------------------------------------------- -# Headless solve (used by run_convergence_dirichlet.py) -# --------------------------------------------------------------------------- - -def solve_dirichlet(elem, mms, L, E, nu, nx, ny, nz, force_field, linear_solver): - """Comme solve_solid, mais avec la baseline = rest_position explicitement - (pas pos0), car ici pos initiale != rest pour les nœuds de bord une fois - l'animation lancée.""" - root = Sofa.Core.Node("root") - dofs, topology = build_dirichlet_scene( - root, mms, elem, L=L, E=E, nu=nu, - nx=nx, ny=ny, nz=nz, with_visual=False, - force_field=force_field, linear_solver=linear_solver - ) - Sofa.Simulation.init(root) - nodes_3d = dofs.rest_position.array().copy() - conn = elem.read_connectivity(topology) - Sofa.Simulation.animate(root, root.dt.value) - pos1 = dofs.position.array().copy() - Sofa.Simulation.unload(root) - 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) - - -# --------------------------------------------------------------------------- -# runSofa entry point (scène visuelle, params.json -> reference block) -# --------------------------------------------------------------------------- - -def createScene(rootNode): - cfg = load_params() - ref = cfg["reference"] - build_dirichlet_scene(rootNode, mms, element_tet, - L=cfg["length"], E=cfg["youngModulus"], - nu=ref["nu"], - nx=ref["nx"], ny=ref["nx"], nz=ref["nx"], - with_visual=True, - force_field=cfg["forceField"], - linear_solver=cfg["linearSolver"]) - return rootNode - - -# --------------------------------------------------------------------------- -# Reference-mesh driver: solve + table + plots (mirrors run_reference_scene) -# --------------------------------------------------------------------------- - -def run_reference_dirichlet(elem, mms): - """Solve au maillage de reference, ecrit la table + les plots PNG - (profil 1D + slices 3D) dans results/, comme run_reference_scene mais - pour le pipeline Dirichlet-only.""" - cfg = load_params() - ref = cfg["reference"] - L, E = cfg["length"], cfg["youngModulus"] - nu = ref["nu"] - nx = ny = nz = ref["nx"] - ff = cfg["forceField"] - ls = cfg["linearSolver"] - - sol = solve_dirichlet(elem, mms, L, E, nu, nx, ny, nz, - force_field=ff, linear_solver=ls) - l2 = elem.compute_l2(sol, mms, L) - h1 = elem.compute_h1(sol, mms, L) - - label = elem.LABEL + " DIRICHLET" - tag = label.replace(" ", "_") - stem = f"{mms.name}_{tag}_nu{nu}_nx{nx}" - - xyz = sol.nodes[:, :3] - write_solution_table(f"solution_{stem}", xyz, - np.column_stack([sol.ux, sol.uy, sol.uz]), - lambda xi, yi, zi: mms.u_ex(xi, yi, zi, L), - RESULTS_DIR, {"L2": l2, "H1_semi": h1}) - plot_solution_profile(f"solution_{stem}", sol, mms, L, nx, ny, nz, - label, nu, l2, h1) - plot_solution_slices(f"fields3D_{stem}", sol, mms, L, nx, ny, nz, - label, nu) - print(f"Dirichlet-only nx={nx} L2={l2:.6e} H1={h1:.6e}") - - -if __name__ == "__main__": - run_reference_dirichlet(element_tet, mms) diff --git a/examples/Freefem/MMS/3D_tetra/sinusoidal_tetra.py b/examples/Freefem/MMS/3D_tetra/sinusoidal_tetra.py deleted file mode 100644 index 984ca8dc..00000000 --- a/examples/Freefem/MMS/3D_tetra/sinusoidal_tetra.py +++ /dev/null @@ -1,13 +0,0 @@ -"""Tet-element entry point for `runSofa`, mirrors sinusoidal.py but loads -the P1 tetrahedron scene (element_tet) instead of the default Q1 hex. - -Usage: runSofa sinusoidal_tet.py -""" - -from sinusoidal import mms -from solid import case_scene, run_reference_scene, element_tet - -createScene = case_scene(mms, element_tet) - -if __name__ == "__main__": - run_reference_scene(element_tet, mms) \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/solid.py b/examples/Freefem/MMS/3D_tetra/solid.py deleted file mode 100644 index 8e5a6abb..00000000 --- a/examples/Freefem/MMS/3D_tetra/solid.py +++ /dev/null @@ -1,316 +0,0 @@ -import json -import numpy as np -import matplotlib.pyplot as plt -import os -import sys - -import Sofa -import Sofa.Core -import Sofa.Simulation -import SofaRuntime - -# 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, element_tet # re-exported for case files -from solid_solution import SolidSolution3D -from output import write_solution_table -from scene import NodalForceAssembler - -RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results") - - -# --------------------------------------------------------------------------- -# Parameters -# --------------------------------------------------------------------------- - -def load_params(path=None): - if path is None: - path = os.path.join(os.path.dirname(os.path.abspath(__file__)), - "params.json") - with open(path) as f: - return json.load(f) - - -# ============== Mesh helpers ============================ - -def get_nodes_3d(L, nx, ny, nz): - dx = L / (nx - 1) - dy = L / (ny - 1) - dz = L / (nz - 1) - pts = [[i*dx, j*dy, k*dz] - for k in range(nz) for j in range(ny) for i in range(nx)] - return np.array(pts, dtype=float) - - -# --------------------------------------------------------------------------- -# SOFA scene -# --------------------------------------------------------------------------- - -def _solid_force_compute(element, mms, L, E, nu, nx, ny, nz): - """Return a `(nodes, topology) -> (N, 3) force array` for NodalForceAssembler.""" - def compute(nodes, topology): - conn = element.read_connectivity(topology) - return element.compute_nodal_forces( - nodes, conn, mms, L, E, nu, nx, ny, nz) - return compute - -def build_solid_scene(rootNode, mms, element, force_field, linear_solver, - L=1.0, E=1e6, nu=0.3, - nx=6, ny=6, nz=6, with_visual=True): - """Build a SOFA scene for `mms` on the 3D `element` strategy. - - Returns (dofs, topology). Nodes and connectivity become available - after `Sofa.Simulation.init(root)` runs, via - `dofs.rest_position.array()` and `element.read_connectivity(topology)`. - - force_field : name of the FEM force field to test - linear_solver : dict {"type": , "parameters": {...}} of the linear solver - """ - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", - "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", - "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.LinearSolver.Iterative", - "Sofa.Component.MechanicalLoad", - "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Grid", - "Sofa.Component.Topology.Container.Dynamic", - "Sofa.Component.Topology.Mapping", - "Sofa.Component.Visual", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_3d = get_nodes_3d(L, nx, ny, nz) - - # Grid must be added *before* Solid: SOFA inits children in insertion - # order, and Solid's topology container resolves `@../Grid/grid.position` - # during its own init — so Grid has to be initialised first. - Grid = rootNode.addChild("Grid") - Grid.addObject("RegularGridTopology", name="grid", - nx=nx, ny=ny, nz=nz, - min=[0.0, 0.0, 0.0], max=[L, L, L]) - - Solid = rootNode.addChild("Solid") - Solid.addObject("StaticSolver", name="staticSolver", printLog=False) - Solid.addObject("NewtonRaphsonSolver", name="newtonSolver", - maxNbIterationsNewton=1, - absoluteResidualStoppingThreshold=1e-10, - printLog=True) - Solid.addObject(linear_solver["type"], name="linearSolver", - printLog=True, - **linear_solver["parameters"]) - - dofs = Solid.addObject("MechanicalObject", name="dofs", template="Vec3d", - 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") - - mms.apply_bcs(Solid, nodes_3d, L) - - # Placeholder force field filled in by the controller after init. - n_nodes = len(nodes_3d) - force_field = Solid.addObject("ConstantForceField", name="MMS_forces", - template="Vec3d", - indices=list(range(n_nodes)), - forces=[[0.0, 0.0, 0.0]] * n_nodes) - - Solid.addObject(NodalForceAssembler( - dofs=dofs, topology=topology, force_field=force_field, - compute_forces=_solid_force_compute(element, mms, L, E, nu, nx, ny, nz), - name="nodalForceAssembler")) - - return dofs, topology - - -# ───────────────────────────────────────────────────────────────────────────── -# Simulation runner -# ───────────────────────────────────────────────────────────────────────────── - -def solve_solid(elem, mms, L, E, nu, nx, ny, nz, force_field, linear_solver): - """Build, init, and run one static step. Returns a SolidSolution3D snapshot.""" - root = Sofa.Core.Node("root") - dofs, topology = build_solid_scene( - root, mms, elem, L=L, E=E, nu=nu, - nx=nx, ny=ny, nz=nz, with_visual=False, - force_field=force_field, linear_solver=linear_solver - ) - Sofa.Simulation.init(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] - return SolidSolution3D(nodes=nodes_3d, conn=conn, ux=ux, uy=uy, uz=uz) - - -# --------------------------------------------------------------------------- -# Output helpers (mirror 2D beam.py) -# --------------------------------------------------------------------------- - -def plot_solution_profile(stem, sol, mms, L, nx, ny, nz, label, nu, l2, h1): - """Save 1-D centerline profiles (ux(x), uy(y), uz(z)) to results/.png.""" - os.makedirs(RESULTS_DIR, exist_ok=True) - xyz = sol.nodes[:, :3] - - mid_i, mid_j, mid_k = nx // 2, ny // 2, nz // 2 - - def nidx(i, j, k): - return i + j * nx + k * nx * ny - - x_fine = np.linspace(0, L, 200) - - line_x = [nidx(i, mid_j, mid_k) for i in range(nx)] - line_y = [nidx(mid_i, j, mid_k) for j in range(ny)] - line_z = [nidx(mid_i, mid_j, k) for k in range(nz)] - - yc, zc = xyz[line_x[0], 1], xyz[line_x[0], 2] - xc2, zc2 = xyz[line_y[0], 0], xyz[line_y[0], 2] - xc3, yc3 = xyz[line_z[0], 0], xyz[line_z[0], 1] - - ux_ex, _, _ = mms.u_ex(x_fine, - np.full_like(x_fine, yc), - np.full_like(x_fine, zc), L) - _, uy_ex, _ = mms.u_ex(np.full_like(x_fine, xc2), - x_fine, - np.full_like(x_fine, zc2), L) - _, _, uz_ex = mms.u_ex(np.full_like(x_fine, xc3), - np.full_like(x_fine, yc3), - x_fine, L) - - fig, axes = plt.subplots(1, 3, figsize=(15, 4)) - for ax, coord, sofa, exact, ylabel, axname in [ - (axes[0], xyz[line_x, 0], sol.ux[line_x], ux_ex, r"$u_x$", "x"), - (axes[1], xyz[line_y, 1], sol.uy[line_y], uy_ex, r"$u_y$", "y"), - (axes[2], xyz[line_z, 2], sol.uz[line_z], uz_ex, r"$u_z$", "z"), - ]: - ax.plot(coord, sofa, "o-", color="tab:green", - label=f"SOFA {label}", ms=5) - ax.plot(x_fine, exact, "--", color="tab:blue", label="MMS exact") - ax.set_xlabel(axname); ax.set_ylabel(ylabel) - ax.legend(); ax.grid(True, alpha=0.3) - fig.suptitle(f"{mms.name} — {label} " - f"nu={nu} nx={nx} |L2={l2:.2e} H1={h1:.2e}") - fig.tight_layout() - fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150) - plt.close(fig) - - -def plot_solution_slices(stem, sol, mms, L, nx, ny, nz, label, nu): - """Save midplane heatmaps (xy at z=mid, xz at y=mid, yz at x=mid) for ux/uy/uz. - - Mirrors the 2-row × 3-column layout of the existing 3D plot_displacement: - each column is a slice plane, each row a displacement component (two per plane). - """ - os.makedirs(RESULTS_DIR, exist_ok=True) - xyz = sol.nodes[:, :3] - mid_i, mid_j, mid_k = nx // 2, ny // 2, nz // 2 - - def nidx(i, j, k): - return i + j * nx + k * nx * ny - - sl_z = np.array([[nidx(i, j, mid_k) for i in range(nx)] for j in range(ny)]) - sl_y = np.array([[nidx(i, mid_j, k) for i in range(nx)] for k in range(nz)]) - sl_x = np.array([[nidx(mid_i, j, k) for j in range(ny)] for k in range(nz)]) - - fig, axes = plt.subplots(2, 3, figsize=(16, 10)) - - # Column 0: z = mid plane — show ux, uy - X, Y = xyz[sl_z, 0], xyz[sl_z, 1] - for row, u, title in [(0, sol.ux, r"$u_x$ z=mid"), - (1, sol.uy, r"$u_y$ z=mid")]: - im = axes[row, 0].pcolormesh(X, Y, u[sl_z], cmap="RdBu_r", shading="auto") - axes[row, 0].set(xlabel="x", ylabel="y", title=title) - plt.colorbar(im, ax=axes[row, 0]) - - # Column 1: y = mid plane — show ux, uz - X, Z = xyz[sl_y, 0], xyz[sl_y, 2] - for row, u, title in [(0, sol.ux, r"$u_x$ y=mid"), - (1, sol.uz, r"$u_z$ y=mid")]: - im = axes[row, 1].pcolormesh(X, Z, u[sl_y], cmap="RdBu_r", shading="auto") - axes[row, 1].set(xlabel="x", ylabel="z", title=title) - plt.colorbar(im, ax=axes[row, 1]) - - # Column 2: x = mid plane — show uy, uz - Y, Z = xyz[sl_x, 1], xyz[sl_x, 2] - for row, u, title in [(0, sol.uy, r"$u_y$ x=mid"), - (1, sol.uz, r"$u_z$ x=mid")]: - im = axes[row, 2].pcolormesh(Y, Z, u[sl_x], cmap="RdBu_r", shading="auto") - axes[row, 2].set(xlabel="y", ylabel="z", title=title) - plt.colorbar(im, ax=axes[row, 2]) - - fig.suptitle(f"Fields 3D — {label} {mms.name} nu={nu} nx={nx}") - fig.tight_layout() - fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150) - plt.close(fig) - - -# ───────────────────────────────────────────────────────────────────────────── -# Single-case driver -# ───────────────────────────────────────────────────────────────────────────── - -def run_reference_scene(elem, mms): - """Solve one MMS case at the reference mesh, write the solution table and plots. - - All parameters come from params.json (top-level + `reference` block). - """ - cfg = load_params() - ref = cfg["reference"] - L, E = cfg["length"], cfg["youngModulus"] - nu = ref["nu"] - nx = ny = nz = ref["nx"] - ff = cfg["forceField"] - ls = cfg["linearSolver"] - - sol = solve_solid(elem, mms, L, E, nu, nx, ny, nz, - force_field=ff, linear_solver=ls) - l2 = elem.compute_l2(sol, mms, L) - h1 = elem.compute_h1(sol, mms, L) - - label = elem.LABEL - tag = label.replace(" ", "_") - stem = f"{mms.name}_{tag}_nu{nu}_nx{nx}" - - xyz = sol.nodes[:, :3] - write_solution_table(f"solution_{stem}", xyz, - np.column_stack([sol.ux, sol.uy, sol.uz]), - lambda xi, yi, zi: mms.u_ex(xi, yi, zi, L), - RESULTS_DIR, {"L2": l2, "H1_semi": h1}) - plot_solution_profile(f"solution_{stem}", sol, mms, L, nx, ny, nz, - label, nu, l2, h1) - plot_solution_slices (f"fields3D_{stem}", sol, mms, L, nx, ny, nz, - label, nu) - - -def case_scene(mms, element): - """Return a `createScene(rootNode)` bound to this MMS and element type. - - All parameters come from params.json (top-level + `reference` block). Each - case file exposes: - createScene = case_scene(mms, element_hex) - so that `runSofa sinusoidal.py` loads the default scene. - """ - def createScene(rootNode): - cfg = load_params() - ref = cfg["reference"] - build_solid_scene(rootNode, mms, element, - L=cfg["length"], E=cfg["youngModulus"], - nu=ref["nu"], - nx=ref["nx"], ny=ref["nx"], nz=ref["nx"], - with_visual=True, force_field=cfg["forceField"], - linear_solver=cfg["linearSolver"]) - return rootNode - return createScene \ No newline at end of file diff --git a/examples/Freefem/MMS/3D_tetra/solid_solution.py b/examples/Freefem/MMS/3D_tetra/solid_solution.py deleted file mode 100644 index 278469f0..00000000 --- a/examples/Freefem/MMS/3D_tetra/solid_solution.py +++ /dev/null @@ -1,14 +0,0 @@ -"""Snapshot of one 3D solve: mesh + connectivity + displacement field.""" - -from dataclasses import dataclass - -import numpy as np - - -@dataclass -class SolidSolution3D: - nodes : np.ndarray # (N, 3) - conn : np.ndarray # (n_hex, 8) hexahedral connectivity from SOFA - ux : np.ndarray # (N,) x-displacement - uy : np.ndarray # (N,) y-displacement - uz : np.ndarray # (N,) z-displacement 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) # ---------------------------------------------------------------------------