From f82726c55b23f41c5d52a42568aa2b170fc62aa1 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 26 May 2026 17:10:21 +0200 Subject: [PATCH 1/2] Add 3D FEM utils to fem.py --- examples/Freefem/MMS/fem.py | 211 ++++++++++++++++++++++++++++++++++++ 1 file changed, 211 insertions(+) diff --git a/examples/Freefem/MMS/fem.py b/examples/Freefem/MMS/fem.py index 146e2678..8288d724 100644 --- a/examples/Freefem/MMS/fem.py +++ b/examples/Freefem/MMS/fem.py @@ -128,6 +128,78 @@ def rule(xe, ye): return rule +# --------------------------------------------------------------------------- +# 3D Q1 reference shape functions on [-1,1]^3. +# Node ordering matches SOFA's RegularGridTopology hex convention: +# 0:(-,-,-) 1:(+,-,-) 2:(+,+,-) 3:(-,+,-) +# 4:(-,-,+) 5:(+,-,+) 6:(+,+,+) 7:(-,+,+) +# --------------------------------------------------------------------------- + +def _shape_hex_q1(xi, eta, zeta): + 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 + + +# --------------------------------------------------------------------------- +# 3D element rules +# +# Same protocol as the 2D rules, lifted one dimension: +# rule(xe, ye, ze) yields (xg, yg, zg, w, N, dN_dx, dN_dy, dN_dz) +# --------------------------------------------------------------------------- + +def hex_q1_rule(n_pts=2): + """Element rule for Q1 hexes: tensor-product Gauss-Legendre over [-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, ye, ze): + 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 @ xe, dN_dxi @ ye, dN_dxi @ ze], + [dN_deta @ xe, dN_deta @ ye, dN_deta @ ze], + [dN_dzeta @ xe, dN_dzeta @ ye, dN_dzeta @ ze], + ]) + detJ = np.linalg.det(J) + Jinv = np.linalg.inv(J) + xg, yg, zg = N @ xe, N @ ye, N @ ze + 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 + yield xg, yg, zg, wi * wj * wk * detJ, N, dN_dx, dN_dy, dN_dz + return rule + + # --------------------------------------------------------------------------- # FEM assembly # --------------------------------------------------------------------------- @@ -211,6 +283,79 @@ def assemble_nodal_forces_2d(f_body, nodes, conn, element_rule): return F +def assemble_traction_3d(traction, nodes, quads, n_pts=2): + """ + Assemble consistent nodal forces from a boundary traction along quad faces: + + F_a += integral_{face} N_a(s,t) * traction(x(s,t), y(s,t), z(s,t)) dA + + Bilinear quad shape functions on each face. The outward unit normal is + baked into `traction` by the caller (one lambda per boundary face) — same + pattern as `assemble_traction_2d`. + + traction : callable (x, y, z) -> (Tx, Ty, Tz) + nodes : (N, 3) array + quads : iterable of (a, b, c, d) node-index 4-tuples on the boundary + n_pts : Gauss points per face direction (default 2 → 2×2 = 4 pts) + """ + if n_pts not in _GAUSS_LEGENDRE_1D: + raise ValueError(f"assemble_traction_3d: {n_pts}-point rule not supported") + xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] + + xyz = np.asarray(nodes)[:, :3] + F = np.zeros((len(xyz), 3)) + for quad in quads: + xe = xyz[quad, 0] + ye = xyz[quad, 1] + ze = xyz[quad, 2] + 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) ]) + # Surface area element from the cross product of the two tangent vectors + t_xi = np.array([dN_dxi @ xe, dN_dxi @ ye, dN_dxi @ ze]) + t_eta = np.array([dN_deta @ xe, dN_deta @ ye, dN_deta @ ze]) + dA = np.linalg.norm(np.cross(t_xi, t_eta)) + xg, yg, zg = N @ xe, N @ ye, N @ ze + Tx, Ty, Tz = traction(xg, yg, zg) + w = wi * wj * dA + for a, node in enumerate(quad): + F[node, 0] += N[a] * Tx * w + F[node, 1] += N[a] * Ty * w + F[node, 2] += N[a] * Tz * w + return F + + +def assemble_nodal_forces_3d(f_body, nodes, conn, element_rule): + """ + Assemble consistent nodal forces on a 3D mesh, agnostic to element type. + + F_a = sum_e integral_{Omega_e} N_a(x,y,z) * f_body(x,y,z) dx dy dz + + f_body : callable (x, y, z) -> (fx, fy, fz) + nodes : (N, 3) array + conn : iterable of node-index lists, one per element + element_rule : rule from hex_q1_rule(n_pts) + """ + xyz = np.asarray(nodes)[:, :3] + F = np.zeros((len(xyz), 3)) + for elem in conn: + xe, ye, ze = xyz[elem, 0], xyz[elem, 1], xyz[elem, 2] + for xg, yg, zg, w, N, _, _, _ in element_rule(xe, ye, ze): + fx, fy, fz = f_body(xg, yg, zg) + for a, node in enumerate(elem): + F[node, 0] += N[a] * fx * w + F[node, 1] += N[a] * fy * w + F[node, 2] += N[a] * fz * w + return F + + # --------------------------------------------------------------------------- # Error norms # --------------------------------------------------------------------------- @@ -296,3 +441,69 @@ def h1_semi_error_2d(nodes, conn, u_h, grad_u_ex, element_rule): (duy_dx_h - G[1, 0])**2 + (duy_dy_h - G[1, 1])**2 ) * w return np.sqrt(total) + + +def l2_error_3d(nodes, conn, u_h, u_ex, element_rule): + """ + Vector L2 error norm on a 3D mesh, element-type agnostic. + + sqrt( integral || u_h(x,y,z) - u_ex(x,y,z) ||^2 dx dy dz ) over the mesh + + nodes : (N, 3) array + conn : iterable of node-index lists, one per element + u_h : (N, 3) array of nodal displacements (ux, uy, uz) + u_ex : callable (x, y, z) -> (ux_ex, uy_ex, uz_ex) + element_rule : rule from hex_q1_rule + """ + xyz = np.asarray(nodes)[:, :3] + u = np.asarray(u_h) + total = 0.0 + for elem in conn: + xe, ye, ze = xyz[elem, 0], xyz[elem, 1], xyz[elem, 2] + u_loc = u[elem] # (n_local, 3) + for xg, yg, zg, w, N, _, _, _ in element_rule(xe, ye, ze): + u_h_g = N @ u_loc # (3,) + ux_ex, uy_ex, uz_ex = u_ex(xg, yg, zg) + ex = u_h_g[0] - ux_ex + ey = u_h_g[1] - uy_ex + ez = u_h_g[2] - uz_ex + total += (ex * ex + ey * ey + ez * ez) * w + return np.sqrt(total) + + +def h1_semi_error_3d(nodes, conn, u_h, grad_u_ex, element_rule): + """ + Vector H1 semi-norm error on a 3D mesh, element-type agnostic. + + sqrt( integral || grad u_h - grad u_ex ||_F^2 dx dy dz ) over the mesh + + nodes : (N, 3) array + conn : iterable of node-index lists, one per element + u_h : (N, 3) array of nodal displacements (ux, uy, uz) + grad_u_ex : callable (x, y, z) -> 3x3 array + [[dux/dx, dux/dy, dux/dz], [duy/dx, ...], [duz/dx, ...]] + element_rule : rule from hex_q1_rule + """ + xyz = np.asarray(nodes)[:, :3] + u = np.asarray(u_h) + total = 0.0 + for elem in conn: + xe, ye, ze = xyz[elem, 0], xyz[elem, 1], xyz[elem, 2] + u_loc = u[elem] # (n_local, 3) + for xg, yg, zg, w, _, dN_dx, dN_dy, dN_dz in element_rule(xe, ye, ze): + dux_dx_h = dN_dx @ u_loc[:, 0] + dux_dy_h = dN_dy @ u_loc[:, 0] + dux_dz_h = dN_dz @ u_loc[:, 0] + duy_dx_h = dN_dx @ u_loc[:, 1] + duy_dy_h = dN_dy @ u_loc[:, 1] + duy_dz_h = dN_dz @ u_loc[:, 1] + duz_dx_h = dN_dx @ u_loc[:, 2] + duz_dy_h = dN_dy @ u_loc[:, 2] + duz_dz_h = dN_dz @ u_loc[:, 2] + G = grad_u_ex(xg, yg, zg) + total += ( + (dux_dx_h - G[0, 0])**2 + (dux_dy_h - G[0, 1])**2 + (dux_dz_h - G[0, 2])**2 + + (duy_dx_h - G[1, 0])**2 + (duy_dy_h - G[1, 1])**2 + (duy_dz_h - G[1, 2])**2 + + (duz_dx_h - G[2, 0])**2 + (duz_dy_h - G[2, 1])**2 + (duz_dz_h - G[2, 2])**2 + ) * w + return np.sqrt(total) From 88afb1a9d910d8eaa19ebaada3283f26a96ccd54 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 26 May 2026 17:12:56 +0200 Subject: [PATCH 2/2] Refactor 3D just like 1D/2D - manufactured_solution.py - solid_solution.py - solid.py - sinus_neumann.py - run_convergence.py - params.json --- .../Freefem/MMS/3D/manufactured_solution.py | 60 ++ examples/Freefem/MMS/3D/params.json | 17 +- examples/Freefem/MMS/3D/run_convergence.py | 137 ++++ examples/Freefem/MMS/3D/sin.py | 71 -- examples/Freefem/MMS/3D/sinus_neumann.py | 113 +++ examples/Freefem/MMS/3D/solid.py | 477 +++++++++++++ examples/Freefem/MMS/3D/solid_solution.py | 14 + examples/Freefem/MMS/3D/utils.py | 659 ------------------ 8 files changed, 810 insertions(+), 738 deletions(-) create mode 100644 examples/Freefem/MMS/3D/manufactured_solution.py create mode 100644 examples/Freefem/MMS/3D/run_convergence.py delete mode 100644 examples/Freefem/MMS/3D/sin.py create mode 100644 examples/Freefem/MMS/3D/sinus_neumann.py create mode 100644 examples/Freefem/MMS/3D/solid.py create mode 100644 examples/Freefem/MMS/3D/solid_solution.py delete mode 100644 examples/Freefem/MMS/3D/utils.py diff --git a/examples/Freefem/MMS/3D/manufactured_solution.py b/examples/Freefem/MMS/3D/manufactured_solution.py new file mode 100644 index 00000000..de39c6d7 --- /dev/null +++ b/examples/Freefem/MMS/3D/manufactured_solution.py @@ -0,0 +1,60 @@ +"""Abstract base class for 3D MMS cases (Cartesian domain [0,L]^3).""" + +from abc import ABC, abstractmethod +import numpy as np + + +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(ABC): + name = None # case identifier (must match the params.json key) + plot_label = None # LaTeX label for the exact solution + + # 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)) + + @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] diff --git a/examples/Freefem/MMS/3D/params.json b/examples/Freefem/MMS/3D/params.json index ce2765ff..6a902e26 100644 --- a/examples/Freefem/MMS/3D/params.json +++ b/examples/Freefem/MMS/3D/params.json @@ -1,13 +1,14 @@ { "length": 1.0, "youngModulus": 1e6, - "RatioPoisson": 0.3, - "nx": 5, - "ny": 5, - "nz": 5, - "nxConvergence": { - "linear_neumann": [5, 15, 26, 36, 46, 56, 66, 76, 86], - "sinus_neumann": [5, 15, 26, 36, 46, 56, 66, 76, 86], - "sinus_dirichlet": [5, 15, 26, 36, 46, 56, 66, 76, 86] + "reference": { + "nx": 6, + "nu": 0.3 + }, + "convergence": { + "nu_values": [0.0, 0.3, 0.49], + "nx_values": { + "sinus_neumann": [4, 6, 8, 12, 16] + } } } diff --git a/examples/Freefem/MMS/3D/run_convergence.py b/examples/Freefem/MMS/3D/run_convergence.py new file mode 100644 index 00000000..84b9dee5 --- /dev/null +++ b/examples/Freefem/MMS/3D/run_convergence.py @@ -0,0 +1,137 @@ +"""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). +""" + +import os +import numpy as np +import matplotlib.pyplot as plt + +from sinus_neumann import mms as sinus_neumann_mms + +from solid import ( + RESULTS_DIR, + load_params, + solve_solid, + element_hex, +) +from plot_utils import annotate_convergence_rates + + +def convergence_study(elem_specs, mms, L, E, nu, nx_values): + """ + Run convergence study 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' + """ + print(f"\n PoissonRatio = {nu}", flush=True) + hdr = (f"{'nx':>5} | {'h':>10} | {'L2':>14} | {'rate_L2':>7} " + f"| {'H1':>14} | {'rate_H1':>7}") + + 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}" + + print(f"\n── {label} {mms.name} nu={nu} ──\n{hdr}", flush=True) + + rows, hs, l2s, h1s = [], [], [], [] + for k, nx in enumerate(nx_values): + ny = nz = nx + h = L / (nx - 1) + sol = solve_solid(elem, mms, L, E, nu, nx, ny, nz) + l2 = elem.compute_l2(sol, mms, L) + h1 = elem.compute_h1(sol, mms, L) + + rate_l2 = (f"{np.log(l2 / l2s[-1]) / np.log(h / hs[-1]):.2f}" + if k > 0 else "") + rate_h1 = (f"{np.log(h1 / h1s[-1]) / np.log(h / hs[-1]):.2f}" + if k > 0 else "") + print(f"{nx:5d} | {h:10.4f} | {l2:14.6e} | {rate_l2:>7} " + f"| {h1:14.6e} | {rate_h1:>7}", flush=True) + rows.append({"nx": nx, "h": h, + "L2": l2, "rate_L2": rate_l2, + "H1": h1, "rate_H1": rate_h1}) + hs.append(h); l2s.append(l2); h1s.append(h1) + + write_convergence_table(stem, rows) + plot_series.append({"label": f"{label} L²", + "errors": l2s, "style": spec["l2_style"]}) + plot_series.append({"label": f"{label} H¹", + "errors": h1s, "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) + + +def write_convergence_table(stem, rows): + """ + Write convergence table to results/.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, ylabel="Error"): + """ + Save log-log convergence plot to results/.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) + + +if __name__ == "__main__": + cfg = load_params() + L = cfg["length"] + E = cfg["youngModulus"] + conv = cfg["convergence"] + + specs = [ + {"elem": element_hex, "label": "Q1 hex", + "l2_style": "bo-", "h1_style": "rs--"}, + ] + + for mms in (sinus_neumann_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) diff --git a/examples/Freefem/MMS/3D/sin.py b/examples/Freefem/MMS/3D/sin.py deleted file mode 100644 index 905288d1..00000000 --- a/examples/Freefem/MMS/3D/sin.py +++ /dev/null @@ -1,71 +0,0 @@ -import numpy as np -from utils import ( - load_params, - run_simulation_3d, - l2_error_3d, - h1_error_3d, - plot_displacement, - plot_convergence, - build_scene_3d, -) - -CASE_NAME = "sinus_neumann" -MODE = "sinus_neumann" - - -def createScene(rootNode): - cfg = load_params() - L = cfg["length"] - E = cfg["youngModulus"] - nu = cfg["RatioPoisson"] - nx = cfg["nx"] - ny = cfg["ny"] - nz = cfg["nz"] - build_scene_3d(rootNode, L=L, E=E, nu=nu, nx=nx, ny=ny, nz=nz, - mode=MODE, visual=True) - return rootNode - - -if __name__ == "__main__": - cfg = load_params() - L = cfg["length"] - E = cfg["youngModulus"] - nu = cfg["RatioPoisson"] - - # ========== ÉTUDE DE CONVERGENCE ========== - nx_list = cfg["nxConvergence"][MODE] - - hs, l2s, h1s = [], [], [] - - for nx in nx_list: - ny = nz = nx - h = L / (nx - 1) - - print(f"\n=== Maillage {nx}×{ny}×{nz} (h={h:.4f}) ===") - - nodes, ux, uy, uz = run_simulation_3d(L, E, nu, nx, ny, nz, mode=MODE) - - l2_err = l2_error_3d(nodes, ux, uy, uz, L, mode=MODE) - h1_err = h1_error_3d(nodes, ux, uy, uz, L, nx, ny, nz, mode=MODE) - - hs.append(h) - l2s.append(l2_err) - h1s.append(h1_err) - - if len(hs) == 1: - print(f" L2 error: {l2_err:.6e}") - print(f" H1 error: {h1_err:.6e}") - else: - l2_rate = np.log(l2s[-2]/l2s[-1]) / np.log(hs[-2]/hs[-1]) - h1_rate = np.log(h1s[-2]/h1s[-1]) / np.log(hs[-2]/hs[-1]) - print(f" L2 error: {l2_err:.6e} (rate {l2_rate:.2f})") - print(f" H1 error: {h1_err:.6e} (rate {h1_rate:.2f})") - - - plot_convergence(hs, l2s, h1s, mode=MODE) - - - nx = nx_list[-1] - ny = nz = nx - nodes, ux, uy, uz = run_simulation_3d(L, E, nu, nx, ny, nz, mode=MODE) - plot_displacement(nodes, ux, uy, uz, L, nx, ny, nz, mode=MODE) diff --git a/examples/Freefem/MMS/3D/sinus_neumann.py b/examples/Freefem/MMS/3D/sinus_neumann.py new file mode 100644 index 00000000..9cfa947e --- /dev/null +++ b/examples/Freefem/MMS/3D/sinus_neumann.py @@ -0,0 +1,113 @@ +""" +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, hex_q1_rule) + + +SINUS_AMPLITUDE = 1e-1 + + +class SinusNeumann(MMSCase3D): + name = "sinus_neumann" + 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)) + + 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) + + # Cauchy momentum: -div sigma = f, with sigma = lam tr(eps) I + 2 mu eps + # for incompressible-style decomposition -((lam+mu) grad(div u) + mu lap u) + 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("sinus_neumann: 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_hex) + + +if __name__ == "__main__": + run_reference_scene(element_hex, mms) diff --git a/examples/Freefem/MMS/3D/solid.py b/examples/Freefem/MMS/3D/solid.py new file mode 100644 index 00000000..ce403450 --- /dev/null +++ b/examples/Freefem/MMS/3D/solid.py @@ -0,0 +1,477 @@ +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 ( + assemble_nodal_forces_3d, + assemble_traction_3d, + l2_error_3d, + h1_semi_error_3d, + hex_q1_rule, +) +from solid_solution import SolidSolution3D + +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) + + +def _boundary_quads(nx, ny, nz): + """Return (xm, xp, ym, yp, zm, zp) quad lists for a structured nx×ny×nz grid. + + Each quad is a 4-tuple of node indices. Node index convention: + `idx(i, j, k) = i + j*nx + k*nx*ny` (matches SOFA RegularGridTopology). + Orientation is consistent within each face but its sign does not affect + `assemble_traction_3d`, which uses |t_xi × t_eta| for the surface area + and takes the outward normal from the caller. + """ + 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 + + +# --------------------------------------------------------------------------- +# Element strategies +# +# Each element type is a thin class declaring: +# LABEL : human-readable identifier (used in plot legends) +# ELEMENT_RULE : an element rule (e.g. hex_q1_rule(2)) +# add_topology(Solid) : wire the SOFA topology and return the topology +# container that holds the connectivity. +# Uses RegularGridTopology under a sibling `Grid` child. +# read_connectivity(topology) : extract the connectivity array from the +# SOFA topology container, post-init. +# +# Assembly / error-norm logic lives once on _ElementBase. Connectivity is +# supplied by the caller (the NodalForceAssembler controller reads it from +# SOFA after init); the element class never holds its own copy. +# --------------------------------------------------------------------------- + +class _ElementBase: + @classmethod + def compute_nodal_forces(cls, nodes_3d, conn, mms, L, E, nu, nx, ny, nz): + xyz = nodes_3d[:, :3] + + F = assemble_nodal_forces_3d( + 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)] + for quads, nrm_x, nrm_y, nrm_z in sides: + F += assemble_traction_3d( + 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) + return F + + @classmethod + def compute_l2(cls, sol, mms, L): + return l2_error_3d( + 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_3d( + 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(_ElementBase): + LABEL = "Q1 hex" + ELEMENT_RULE = staticmethod(hex_q1_rule(2)) # used for L²/H¹ error norms + + @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() + + +# --------------------------------------------------------------------------- +# Force assembly controller +# +# Same pattern as 2D: defer assembly to onSimulationInitDoneEvent so that +# SOFA topology components (RegularGridTopology, +# HexahedronSetTopologyContainer, ...) have been initialised and connectivity +# can be read off the topology container. +# --------------------------------------------------------------------------- + +class NodalForceAssembler(Sofa.Core.Controller): + def __init__(self, dofs, topology, force_field, element, mms, + L, E, nu, nx, ny, nz, *args, **kwargs): + super().__init__(*args, **kwargs) + self.dofs = dofs + self.topology = topology + self.force_field = force_field + self.element = element + self.mms = mms + self.L, self.E, self.nu = L, E, nu + self.nx, self.ny, self.nz = nx, ny, nz + + def onSimulationInitDoneEvent(self, event): + nodes = self.dofs.rest_position.array().copy() + conn = self.element.read_connectivity(self.topology) + F = self.element.compute_nodal_forces( + nodes, conn, self.mms, + self.L, self.E, self.nu, self.nx, self.ny, self.nz) + with self.force_field.forces.writeableArray() as forces: + forces[:] = F + + +# --------------------------------------------------------------------------- +# SOFA scene +# --------------------------------------------------------------------------- + +def build_solid_scene(rootNode, mms, element, 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)`. + """ + rootNode.addObject("RequiredPlugin", pluginName=[ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.LinearSolver.Direct", + "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=False) + Solid.addObject("SparseLDLSolver", name="linearSolver", + template="CompressedRowSparseMatrixd") + + 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("LinearSmallStrainFEMForceField", 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, + element=element, mms=mms, + L=L, E=E, nu=nu, nx=nx, ny=ny, nz=nz, + name="nodalForceAssembler")) + + return dofs, topology + + +# ───────────────────────────────────────────────────────────────────────────── +# Simulation runner +# ───────────────────────────────────────────────────────────────────────────── + +def solve_solid(elem, mms, L, E, nu, nx, ny, nz): + """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 + ) + 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 write_solution_table(stem, x, y, z, ux_h, uy_h, uz_h, u_ex, error_dict): + """ + Write per-node solution table and error summary to results/.txt. + + u_ex : callable (x, y, z) -> (ux_ex, uy_ex, uz_ex) + error_dict : ordered dict of {label: value} for error summary lines + """ + os.makedirs(RESULTS_DIR, exist_ok=True) + path = os.path.join(RESULTS_DIR, f"{stem}.txt") + with open(path, "w") as f: + f.write(f"{'x':>10} | {'y':>10} | {'z':>10} | " + f"{'ux_h':>15} | {'uy_h':>15} | {'uz_h':>15} | " + f"{'ux_ex':>15} | {'uy_ex':>15} | {'uz_ex':>15} | " + f"{'err_x':>15} | {'err_y':>15} | {'err_z':>15}\n") + f.write("-" * 184 + "\n") + for xi, yi, zi, uxi, uyi, uzi in zip(x, y, z, ux_h, uy_h, uz_h): + uxe, uye, uze = u_ex(xi, yi, zi) + f.write(f"{xi:10.4f} | {yi:10.4f} | {zi:10.4f} | " + f"{uxi:15.6e} | {uyi:15.6e} | {uzi:15.6e} | " + f"{uxe:15.6e} | {uye:15.6e} | {uze:15.6e} | " + f"{abs(uxi - uxe):15.6e} | {abs(uyi - uye):15.6e} | " + f"{abs(uzi - uze):15.6e}\n") + f.write("\n") + for label, val in error_dict.items(): + f.write(f"{label:12s} = {val:.6e}\n") + + +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"] + + sol = solve_solid(elem, mms, L, E, nu, nx, ny, nz) + 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[:, 0], xyz[:, 1], xyz[:, 2], + sol.ux, sol.uy, sol.uz, + lambda xi, yi, zi: mms.u_ex(xi, yi, zi, L), + {"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 sinus_neumann.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) + return rootNode + return createScene + + +# ───────────────────────────────────────────────────────────────────────────── +# Element instances +# ───────────────────────────────────────────────────────────────────────────── + +element_hex = _HexElement() diff --git a/examples/Freefem/MMS/3D/solid_solution.py b/examples/Freefem/MMS/3D/solid_solution.py new file mode 100644 index 00000000..278469f0 --- /dev/null +++ b/examples/Freefem/MMS/3D/solid_solution.py @@ -0,0 +1,14 @@ +"""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/3D/utils.py b/examples/Freefem/MMS/3D/utils.py deleted file mode 100644 index 3b926092..00000000 --- a/examples/Freefem/MMS/3D/utils.py +++ /dev/null @@ -1,659 +0,0 @@ -import json -import os -import numpy as np -import matplotlib.pyplot as plt -import Sofa -import Sofa.Core -import Sofa.Simulation - -RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results_3d") -SINUS_AMPLITUDE = 1e-1 - - -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) - - -def lame(E, nu): - lam = E * nu / ((1 + nu) * (1 - 2 * nu)) - mu = E / (2 * (1 + nu)) - return lam, mu - - -def node_idx(i, j, k, nx, ny): - return i + j * nx + k * nx * ny - - -def _build_element_connectivity(nx, ny, nz): - i = np.arange(nx-1); j = np.arange(ny-1); k = np.arange(nz-1) - ii, jj, kk = np.meshgrid(i, j, k, indexing='ij') - ii, jj, kk = ii.ravel(), jj.ravel(), kk.ravel() - return np.column_stack([ - ii + jj*nx + kk*nx*ny, - ii+1 + jj*nx + kk*nx*ny, - ii+1 + (jj+1)*nx + kk*nx*ny, - ii + (jj+1)*nx + kk*nx*ny, - ii + jj*nx + (kk+1)*nx*ny, - ii+1 + jj*nx + (kk+1)*nx*ny, - ii+1 + (jj+1)*nx + (kk+1)*nx*ny, - ii + (jj+1)*nx + (kk+1)*nx*ny, - ]) # (n_e, 8) - - -def u_ex_linear(x, y, z, L): - return x / L, y / L, z / L - - -def stress_linear(E, nu, L): - lam, mu = lame(E, nu) - eps = np.diag([1/L, 1/L, 1/L]) - return lam * np.trace(eps) * np.eye(3) + 2 * mu * eps - - -def traction_linear(x, y, z, nx_c, ny_c, nz_c, E, nu, L): - t = stress_linear(E, nu, L) @ np.array([nx_c, ny_c, nz_c]) - return t[0], t[1], t[2] - - - -def u_ex_sinus_neumann(x, y, z, L): - return (SINUS_AMPLITUDE * np.sin(np.pi*x/L)*np.sin(np.pi*y/L), - SINUS_AMPLITUDE * np.sin(np.pi*y/L)*np.sin(np.pi*z/L), - SINUS_AMPLITUDE * np.sin(np.pi*z/L)*np.sin(np.pi*x/L)) - -def _sn_dux_dx(x,y,z,L): return (np.pi/L)*np.cos(np.pi*x/L)*np.sin(np.pi*y/L) -def _sn_dux_dy(x,y,z,L): return (np.pi/L)*np.sin(np.pi*x/L)*np.cos(np.pi*y/L) -def _sn_dux_dz(x,y,z,L): return np.zeros_like(x) if hasattr(x,'__len__') else 0.0 -def _sn_duy_dx(x,y,z,L): return np.zeros_like(x) if hasattr(x,'__len__') else 0.0 -def _sn_duy_dy(x,y,z,L): return (np.pi/L)*np.cos(np.pi*y/L)*np.sin(np.pi*z/L) -def _sn_duy_dz(x,y,z,L): return (np.pi/L)*np.sin(np.pi*y/L)*np.cos(np.pi*z/L) -def _sn_duz_dx(x,y,z,L): return (np.pi/L)*np.sin(np.pi*z/L)*np.cos(np.pi*x/L) -def _sn_duz_dy(x,y,z,L): return np.zeros_like(x) if hasattr(x,'__len__') else 0.0 -def _sn_duz_dz(x,y,z,L): return (np.pi/L)*np.cos(np.pi*z/L)*np.sin(np.pi*x/L) - -def stress_sinus(x, y, z, E, nu, L): - lam, mu = lame(E, nu) - exx = _sn_dux_dx(x,y,z,L); eyy = _sn_duy_dy(x,y,z,L); ezz = _sn_duz_dz(x,y,z,L) - exy = 0.5*(_sn_dux_dy(x,y,z,L) + _sn_duy_dx(x,y,z,L)) - exz = 0.5*(_sn_dux_dz(x,y,z,L) + _sn_duz_dx(x,y,z,L)) - eyz = 0.5*(_sn_duy_dz(x,y,z,L) + _sn_duz_dy(x,y,z,L)) - tr = exx + eyy + ezz - return SINUS_AMPLITUDE * np.array([ - [lam*tr+2*mu*exx, 2*mu*exy, 2*mu*exz ], - [2*mu*exy, lam*tr+2*mu*eyy, 2*mu*eyz ], - [2*mu*exz, 2*mu*eyz, lam*tr+2*mu*ezz], - ]) - -def traction_sinus(x, y, z, nx_c, ny_c, nz_c, E, nu, L): - t = stress_sinus(x, y, z, E, nu, L) @ np.array([nx_c, ny_c, nz_c]) - return t[0], t[1], t[2] - -def body_force_sinus(x, y, z, E, nu, L): - lam, mu = lame(E, nu) - 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) - lap_ux = -2*p**2 * sx*sy - lap_uy = -2*p**2 * sy*sz - lap_uz = -2*p**2 * sz*sx - 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 = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dx - mu*lap_ux) - fy = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dy - mu*lap_uy) - fz = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dz - mu*lap_uz) - return fx, fy, fz - - - -def u_ex_sinus_dirichlet(x, y, z, L): - val = SINUS_AMPLITUDE * np.sin(np.pi*x/L)*np.sin(np.pi*y/L)*np.sin(np.pi*z/L) - return val, val, val - -def _sd_dx(x,y,z,L): return (np.pi/L)*np.cos(np.pi*x/L)*np.sin(np.pi*y/L)*np.sin(np.pi*z/L) -def _sd_dy(x,y,z,L): return (np.pi/L)*np.sin(np.pi*x/L)*np.cos(np.pi*y/L)*np.sin(np.pi*z/L) -def _sd_dz(x,y,z,L): return (np.pi/L)*np.sin(np.pi*x/L)*np.sin(np.pi*y/L)*np.cos(np.pi*z/L) - -def body_force_sinus_dirichlet(x, y, z, E, nu, L): - lam, mu = lame(E, nu) - 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) - - lap_s3 = -3*p**2 * sx*sy*sz - d_divu_dx = p**2*(-sx*sy*sz + cx*cy*sz + cx*sy*cz) - d_divu_dy = p**2*( cx*cy*sz - sx*sy*sz + sx*cy*cz) - d_divu_dz = p**2*( cx*sy*cz + sx*cy*cz - sx*sy*sz) - fx = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dx - mu*lap_s3) - fy = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dy - mu*lap_s3) - fz = SINUS_AMPLITUDE * (-(lam+mu)*d_divu_dz - mu*lap_s3) - return fx, fy, fz - -def get_u_ex(mode): - return {'linear_neumann': u_ex_linear, - 'sinus_neumann': u_ex_sinus_neumann, - 'sinus_dirichlet': u_ex_sinus_dirichlet}[mode] - -def get_derivatives(mode): - if mode == 'linear_neumann': - one = lambda x,y,z,L: np.ones_like(x)/L if hasattr(x,'__len__') else 1/L - zer = lambda x,y,z,L: np.zeros_like(x) if hasattr(x,'__len__') else 0.0 - return {'dux_dx':one,'dux_dy':zer,'dux_dz':zer, - 'duy_dx':zer,'duy_dy':one,'duy_dz':zer, - 'duz_dx':zer,'duz_dy':zer,'duz_dz':one} - elif mode == 'sinus_neumann': - def _s(fn): return lambda x,y,z,L,_f=fn: _f(x,y,z,L)*SINUS_AMPLITUDE - return {'dux_dx':_s(_sn_dux_dx),'dux_dy':_s(_sn_dux_dy),'dux_dz':_s(_sn_dux_dz), - 'duy_dx':_s(_sn_duy_dx),'duy_dy':_s(_sn_duy_dy),'duy_dz':_s(_sn_duy_dz), - 'duz_dx':_s(_sn_duz_dx),'duz_dy':_s(_sn_duz_dy),'duz_dz':_s(_sn_duz_dz)} - else: # sinus_dirichlet : ux=uy=uz=sin*sin*sin - def _s(fn): return lambda x,y,z,L,_f=fn: _f(x,y,z,L)*SINUS_AMPLITUDE - sd = _s(_sd_dx); se = _s(_sd_dy); sf = _s(_sd_dz) - return {'dux_dx':sd,'dux_dy':se,'dux_dz':sf, - 'duy_dx':sd,'duy_dy':se,'duy_dz':sf, - 'duz_dx':sd,'duz_dy':se,'duz_dz':sf} - - - -def _hexa_shape(xi, eta, zeta): - 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 - - -# Precomputed shape functions and derivatives at all 8 Gauss points (2-point rule, gw=1) -_GP = np.array([-1/np.sqrt(3), 1/np.sqrt(3)]) -_N_GP = np.empty((8, 8)) # (n_gp, n_nodes) -_dN_GP = np.empty((8, 3, 8)) # (n_gp, 3, n_nodes) -for _g, (_xi, _eta, _zeta) in enumerate( - [(_xi, _eta, _zeta) for _xi in _GP for _eta in _GP for _zeta in _GP]): - _N, _dNxi, _dNeta, _dNzeta = _hexa_shape(_xi, _eta, _zeta) - _N_GP[_g] = _N - _dN_GP[_g] = np.array([_dNxi, _dNeta, _dNzeta]) - - -def _batch_det3(J): - """Analytical determinant for a batch of (n, 3, 3) matrices.""" - return (J[:,0,0] * (J[:,1,1]*J[:,2,2] - J[:,1,2]*J[:,2,1]) - - J[:,0,1] * (J[:,1,0]*J[:,2,2] - J[:,1,2]*J[:,2,0]) - + J[:,0,2] * (J[:,1,0]*J[:,2,1] - J[:,1,1]*J[:,2,0])) - - -def _batch_detinv3(J): - """Analytical det and inverse for (n, 3, 3) matrices via cofactors.""" - a=J[:,0,0]; b=J[:,0,1]; c=J[:,0,2] - d=J[:,1,0]; e=J[:,1,1]; f=J[:,1,2] - g=J[:,2,0]; h=J[:,2,1]; k=J[:,2,2] - c00=e*k-f*h; c10=f*g-d*k; c20=d*h-e*g - c01=c*h-b*k; c11=a*k-c*g; c21=b*g-a*h - c02=b*f-c*e; c12=c*d-a*f; c22=a*e-b*d - det = a*c00 + b*c10 + c*c20 - s = 1.0 / det - Jinv = np.empty_like(J) - Jinv[:,0,0]=c00*s; Jinv[:,0,1]=c01*s; Jinv[:,0,2]=c02*s - Jinv[:,1,0]=c10*s; Jinv[:,1,1]=c11*s; Jinv[:,1,2]=c12*s - Jinv[:,2,0]=c20*s; Jinv[:,2,1]=c21*s; Jinv[:,2,2]=c22*s - return det, Jinv - - -def compute_body_force_nodal(nodes, nx, ny, nz, E, nu, L, body_force_fn): - all_idx = _build_element_connectivity(nx, ny, nz) # (n_e, 8) - xe_all = nodes[all_idx] # (n_e, 8, 3) - F = np.zeros((len(nodes), 3)) - - for g in range(8): - J_all = _dN_GP[g] @ xe_all # (n_e, 3, 3) - detJ_all = _batch_det3(J_all) # (n_e,) - xg_all = xe_all.transpose(0, 2, 1) @ _N_GP[g] # (n_e, 3) - - fx, fy, fz = body_force_fn(xg_all[:,0], xg_all[:,1], xg_all[:,2], E, nu, L) - f_all = np.stack([fx, fy, fz], axis=1) # (n_e, 3) - w_all = np.abs(detJ_all) # gw=1 for all points - - # contrib[e, local_node, d] = N_g[local_node] * f_all[e,d] * w_all[e] - contrib = (_N_GP[g][np.newaxis, :, np.newaxis] - * f_all[:, np.newaxis, :] - * w_all[:, np.newaxis, np.newaxis]) # (n_e, 8, 3) - for d in range(3): - F[:, d] += np.bincount(all_idx.ravel(), - weights=contrib[:, :, d].ravel(), - minlength=len(nodes)) - return F - - - -def compute_surface_traction_nodal(nodes, quads, E, nu, L, traction_fn): - - gp = np.array([-1/np.sqrt(3), 1/np.sqrt(3)]) - gw = np.array([1.0, 1.0]) - F = np.zeros((len(nodes), 3)) - centroid = nodes.mean(axis=0) - - for quad in quads: - x = nodes[quad] - dN_dxi0 = np.array([-1., 1., 1., -1.]) / 4 - dN_deta0 = np.array([-1., -1., 1., 1.]) / 4 - J0 = np.cross(dN_dxi0 @ x, dN_deta0 @ x) - sign = 1.0 if np.dot(J0, x.mean(axis=0) - centroid) > 0 else -1.0 - - for gi, xi in enumerate(gp): - for gj, eta in enumerate(gp): - N = np.array([(1-xi)*(1-eta),(1+xi)*(1-eta), - (1+xi)*(1+eta),(1-xi)*(1+eta)]) / 4 - dN_dxi = np.array([-(1-eta), (1-eta), (1+eta),-(1+eta)]) / 4 - dN_deta = np.array([-(1-xi), -(1+xi), (1+xi), (1-xi) ]) / 4 - J_vec = np.cross(dN_dxi @ x, dN_deta @ x) - Jmag = np.linalg.norm(J_vec) - n_hat = sign * J_vec / Jmag - xg = N @ x - T = np.array(traction_fn(*xg, *n_hat, E, nu, L)) - F[quad] += np.outer(N, T) * (gw[gi]*gw[gj]*Jmag) - - residual = np.abs(F.sum(axis=0)) - return F - - -class MMSForceController(Sofa.Core.Controller): - def __init__(self, dofs, cff, quad_topo, nx, ny, nz, - E, nu, L, mode, *args, **kwargs): - super().__init__(*args, **kwargs) - self.dofs = dofs - self.cff = cff - self.quad_topo = quad_topo - self.nx, self.ny, self.nz = nx, ny, nz - self.E, self.nu, self.L = E, nu, L - self.mode = mode - - def onSimulationInitDoneEvent(self, event): - nodes = self.dofs.position.array().copy() - quads = np.array(self.quad_topo.quads.array()) - E, nu, L = self.E, self.nu, self.L - nx, ny, nz = self.nx, self.ny, self.nz - mode = self.mode - - if mode == 'linear_neumann': - F = compute_surface_traction_nodal( - nodes, quads, E, nu, L, traction_linear) - elif mode == 'sinus_neumann': - F_body = compute_body_force_nodal(nodes, nx, ny, nz, E, nu, L, body_force_sinus) - F_surf = compute_surface_traction_nodal(nodes, quads, E, nu, L, traction_sinus) - print(f" [body] résidu global F_body : {np.abs(F_body.sum(axis=0))}") - print(f" [total] résidu F_body+F_surf : {np.abs((F_body + F_surf).sum(axis=0))}") - F = F_body + F_surf - elif mode == 'sinus_dirichlet': - F = compute_body_force_nodal( - nodes, nx, ny, nz, E, nu, L, body_force_sinus_dirichlet) - else: - raise ValueError(f"Unknown: {mode}") - - self.cff.forces.value = F.tolist() - - -# =========================== SOFA scene ================================================ - -def build_scene_3d(rootNode, L=1.0, E=1e6, nu=0.3, nx=6, ny=6, nz=6, - mode='linear_neumann', visual=False): - plugins = [ - "Elasticity", - "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", - "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", - "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Dynamic", - "Sofa.Component.Topology.Container.Grid", - "Sofa.Component.Topology.Mapping", - "Sofa.Component.LinearSolver.Iterative", - "Sofa.Component.LinearSolver.Preconditioner", - "Sofa.Component.LinearSystem", - "Sofa.Component.Mapping.Linear" - ] - if visual: - plugins += ["Sofa.Component.Visual","Sofa.GL.Component.Rendering3D"] - - rootNode.addObject('RequiredPlugin', pluginName=plugins) - rootNode.addObject('DefaultAnimationLoop') - if visual: - rootNode.addObject('DefaultVisualManagerLoop') - rootNode.addObject('VisualStyle', - displayFlags='showVisualModels showWireframe') - rootNode.gravity.value = [0,0,0] - rootNode.dt.value = 1.0 - - regGrid = rootNode.addChild('GridTopology') - regGrid.addObject('RegularGridTopology', name="grid", - nx=nx, ny=ny, nz=nz, - min=[0.,0.,0.], max=[L,L,L]) - - solid = rootNode.addChild('Solid3D') - solid.addObject('NewtonRaphsonSolver', name="newtonSolver", - printLog=False, warnWhenLineSearchFails=True - , maxNbIterationsNewton=2 - , relativeSuccessiveStoppingThreshold=0 - , absoluteResidualStoppingThreshold=0 - , absoluteEstimateDifferenceThreshold=0 - , relativeInitialStoppingThreshold=0 - , relativeEstimateDifferenceThreshold=0) - solid.addObject("MatrixLinearSystem", name="ssorSystem", template="CompressedRowSparseMatrixd") - solid.addObject("SSORPreconditioner", name="preconditioner", printLog=False) - solid.addObject("PreconditionedMatrixFreeSystem", name="solverSystem", template="GraphScattered", preconditionerSystem="@ssorSystem") - solid.addObject("PCGLinearSolver", name="linearSolver", template="GraphScattered" - , iterations=3000, tolerance=1e-15, preconditioner="@preconditioner", printLog=True) - #solid.addObject('SparseLDLSolver', name="linearSolver", template="CompressedRowSparseMatrixd") - solid.addObject('StaticSolver', name="staticSolver", - newtonSolver="@newtonSolver", linearSolver="@linearSolver") - solid.addObject('HexahedronSetTopologyContainer', name="topology", - src="@../GridTopology/grid") - solid.addObject('HexahedronSetTopologyModifier') - dofs = solid.addObject('MechanicalObject', name="dofs", - template="Vec3d", src="@topology") - solid.addObject('LinearSmallStrainFEMForceField', name="FEM", - template="Vec3d", youngModulus=E, poissonRatio=nu, - topology="@topology") - - eps = 1e-8 - if mode == 'sinus_dirichlet': - # u_ex=0 - for name_bc, box in [ - ("xm", [-eps, -eps, -eps, eps, L+eps, L+eps]), - ("xp", [L-eps, -eps, -eps, L+eps, L+eps, L+eps]), - ("ym", [-eps, -eps, -eps, L+eps, eps, L+eps]), - ("yp", [-eps, L-eps, -eps, L+eps, L+eps, L+eps]), - ("zm", [-eps, -eps, -eps, L+eps, L+eps, eps ]), - ("zp", [-eps, -eps, L-eps, L+eps, L+eps, L+eps]), - ]: - solid.addObject('BoxROI', name=f"roi_{name_bc}", box=box) - solid.addObject('FixedProjectiveConstraint', - name=f"bc_{name_bc}", - indices=f"@roi_{name_bc}.indices") - else: - solid.addObject('BoxROI', name="fixA", - box=[-eps,-eps,-eps, eps,eps,eps]) - solid.addObject('FixedProjectiveConstraint', name="bcA", - indices="@fixA.indices") - solid.addObject('BoxROI', name="fixB", - box=[L-eps,-eps,-eps, L+eps,eps,eps]) - solid.addObject('PartialFixedProjectiveConstraint', name="bcB", - indices="@fixB.indices", fixedDirections=[0,1,1]) - solid.addObject('BoxROI', name="fixC", - box=[-eps,L-eps,-eps, eps,L+eps,eps]) - solid.addObject('PartialFixedProjectiveConstraint', name="bcC", - indices="@fixC.indices", fixedDirections=[0,0,1]) - - # ======================== Forces MMS ========================================== - n_nodes = nx * ny * nz - cff = solid.addObject('ConstantForceField', name="MMS_forces", - template="Vec3d", - indices=list(range(n_nodes)), - forces=[[0.,0.,0.]] * n_nodes) - - surf = solid.addChild('Surface') - quad_topo = surf.addObject('QuadSetTopologyContainer', name="surfTopo") - surf.addObject('QuadSetTopologyModifier') - surf.addObject('Hexa2QuadTopologicalMapping', - input="@../topology", output="@surfTopo") - if visual: - surf.addObject('OglModel', name="ogl", src="@surfTopo", - color=[0.2,0.6,1.0,0.9]) - surf.addObject('IdentityMapping') - - solid.addObject(MMSForceController( - dofs, cff, quad_topo, nx, ny, nz, E, nu, L, mode, - name="MMSForceController")) - - return dofs - - -def run_simulation_3d(L, E, nu, nx, ny, nz, mode='linear_neumann'): - root = Sofa.Core.Node("root") - root.dt.value = 1.0 - dofs = build_scene_3d(root, L=L, E=E, nu=nu, - nx=nx, ny=ny, nz=nz, mode=mode, visual=False) - Sofa.Simulation.init(root) - pos_init = dofs.position.array().copy() - Sofa.Simulation.animate(root, root.dt.value) - pos_final = dofs.position.array().copy() - Sofa.Simulation.unload(root) - ux = pos_final[:,0] - pos_init[:,0] - uy = pos_final[:,1] - pos_init[:,1] - uz = pos_final[:,2] - pos_init[:,2] - return pos_init, ux, uy, uz - - - -def l2_error_3d(nodes, ux, uy, uz, L, mode): - u_ex = get_u_ex(mode) - ex_x, ex_y, ex_z = u_ex(nodes[:,0], nodes[:,1], nodes[:,2], L) - return np.sqrt(np.mean((ux-ex_x)**2 + (uy-ex_y)**2 + (uz-ex_z)**2)) - - -def h1_error_3d(nodes, ux, uy, uz, L, nx, ny, nz, mode): - x, y, z = nodes[:,0], nodes[:,1], nodes[:,2] - u_ex_fn = get_u_ex(mode) - ex_x, ex_y, ex_z = u_ex_fn(x, y, z, L) - - h = L / (nx - 1) - err_u2 = np.sum((ux-ex_x)**2 + (uy-ex_y)**2 + (uz-ex_z)**2) * h**3 - - d = get_derivatives(mode) - all_idx = _build_element_connectivity(nx, ny, nz) # (n_e, 8) - xe_all = nodes[all_idx] # (n_e, 8, 3) - ue_all = np.column_stack([ux, uy, uz])[all_idx] # (n_e, 8, 3) - - err_grad2 = 0.0 - for g in range(8): - J_all = _dN_GP[g] @ xe_all # (n_e, 3, 3) - detJ_all, Jinv_all = _batch_detinv3(J_all) # (n_e,), (n_e, 3, 3) - - # dN_phys_all[e,i,j] = ∂N_j/∂x_i - dN_phys_all = Jinv_all @ _dN_GP[g] # (n_e, 3, 8) - # grad_uh_all[e,i,j] = ∂u_j/∂x_i - grad_uh_all = dN_phys_all @ ue_all # (n_e, 3, 3) - - xg_all = xe_all.transpose(0, 2, 1) @ _N_GP[g] # (n_e, 3) - xg0, xg1, xg2 = xg_all[:,0], xg_all[:,1], xg_all[:,2] - - # build (3,3,n_e) then transpose to (n_e,3,3) - grad_uex_all = np.array([ - [d['dux_dx'](xg0,xg1,xg2,L), d['duy_dx'](xg0,xg1,xg2,L), d['duz_dx'](xg0,xg1,xg2,L)], - [d['dux_dy'](xg0,xg1,xg2,L), d['duy_dy'](xg0,xg1,xg2,L), d['duz_dy'](xg0,xg1,xg2,L)], - [d['dux_dz'](xg0,xg1,xg2,L), d['duy_dz'](xg0,xg1,xg2,L), d['duz_dz'](xg0,xg1,xg2,L)], - ]).transpose(2, 0, 1) # (n_e, 3, 3) - - err_e = np.sum((grad_uh_all - grad_uex_all)**2, axis=(1, 2)) # (n_e,) - err_grad2 += np.dot(err_e, np.abs(detJ_all)) # gw=1 - - return np.sqrt(err_u2 + err_grad2) - -def plot_convergence(hs, l2s, h1s, mode, out_dir=None): - if out_dir is None: - out_dir = RESULTS_DIR - os.makedirs(out_dir, exist_ok=True) - - hs = np.array(hs); l2s = np.array(l2s); h1s = np.array(h1s) - ref2 = l2s[-1] * (hs / hs[-1])**2 - ref1 = h1s[-1] * (hs / hs[-1])**1 - - fig, axes = plt.subplots(1, 2, figsize=(12, 5)) - - axes[0].loglog(hs, l2s, 'bo-', label='L2 SOFA', lw=2, markersize=7) - axes[0].loglog(hs, ref2, 'r--', label='ref $h^2$', lw=1.5) - axes[0].set(xlabel='h', ylabel='erreur L2', - title=f'Convergence L2 [{mode}]') - axes[0].legend(); axes[0].grid(True, which='both', alpha=0.3) - if len(hs) > 1: - for k, o in enumerate(np.diff(np.log(l2s))/np.diff(np.log(hs))): - axes[0].annotate(f'{o:.2f}', - xy=(np.sqrt(hs[k]*hs[k+1]), np.sqrt(l2s[k]*l2s[k+1])), - fontsize=9, color='blue') - - axes[1].loglog(hs, h1s, 'gs-', label='H1 SOFA', lw=2, markersize=7) - axes[1].loglog(hs, ref1, 'r--', label='ref $h^1$', lw=1.5) - axes[1].set(xlabel='h', ylabel='erreur H1', - title=f'Convergence H1 [{mode}]') - axes[1].legend(); axes[1].grid(True, which='both', alpha=0.3) - if len(hs) > 1: - for k, o in enumerate(np.diff(np.log(h1s))/np.diff(np.log(hs))): - axes[1].annotate(f'{o:.2f}', - xy=(np.sqrt(hs[k]*hs[k+1]), np.sqrt(h1s[k]*h1s[k+1])), - fontsize=9, color='green') - - plt.suptitle(f'MMS 3D Hexa Q1 — {mode}', fontsize=13) - plt.tight_layout() - out = os.path.join(out_dir, f'convergence_{mode}.png') - plt.savefig(out, dpi=150); plt.close() - - - -def plot_displacement(nodes, ux, uy, uz, L, nx, ny, nz, mode, out_dir=None): - if out_dir is None: - out_dir = RESULTS_DIR - os.makedirs(out_dir, exist_ok=True) - - def nidx(i,j,k): return node_idx(i,j,k,nx,ny) - mid_i=nx//2; mid_j=ny//2; mid_k=nz//2 - x_fine = np.linspace(0, L, 200) - u_ex_fn = get_u_ex(mode) - - fig, axes = plt.subplots(1, 3, figsize=(15,5)) - - line_x = [nidx(i,mid_j,mid_k) for i in range(nx)] - y_mid=nodes[line_x[0],1]; z_mid=nodes[line_x[0],2] - ux_ex,_,_ = u_ex_fn(x_fine, np.full_like(x_fine,y_mid), np.full_like(x_fine,z_mid), L) - axes[0].plot(nodes[line_x,0], ux[line_x], 'bo-', label='SOFA', markersize=7) - axes[0].plot(x_fine, ux_ex, 'r--', label='MMS exact', lw=2) - axes[0].set(xlabel='x', ylabel='$u_x$', - title=f'$u_x$(x) y={y_mid:.2f}, z={z_mid:.2f}') - axes[0].legend(); axes[0].grid(alpha=0.3) - - line_y = [nidx(mid_i,j,mid_k) for j in range(ny)] - x_mid2=nodes[line_y[0],0]; z_mid2=nodes[line_y[0],2] - _,uy_ex,_ = u_ex_fn(np.full_like(x_fine,x_mid2), x_fine, np.full_like(x_fine,z_mid2), L) - axes[1].plot(nodes[line_y,1], uy[line_y], 'go-', label='SOFA', markersize=7) - axes[1].plot(x_fine, uy_ex, 'r--', label='MMS exact', lw=2) - axes[1].set(xlabel='y', ylabel='$u_y$', - title=f'$u_y$(y) x={x_mid2:.2f}, z={z_mid2:.2f}') - axes[1].legend(); axes[1].grid(alpha=0.3) - - line_z = [nidx(mid_i,mid_j,k) for k in range(nz)] - x_mid3=nodes[line_z[0],0]; y_mid3=nodes[line_z[0],1] - _,_,uz_ex = u_ex_fn(np.full_like(x_fine,x_mid3), np.full_like(x_fine,y_mid3), x_fine, L) - axes[2].plot(nodes[line_z,2], uz[line_z], 'ms-', label='SOFA', markersize=7) - axes[2].plot(x_fine, uz_ex, 'r--', label='MMS exact', lw=2) - axes[2].set(xlabel='z', ylabel='$u_z$', - title=f'$u_z$(z) x={x_mid3:.2f}, y={y_mid3:.2f}') - axes[2].legend(); axes[2].grid(alpha=0.3) - - plt.suptitle(f'MMS 3D Hexa Q1 [{mode}] — 1D ({nx}×{ny}×{nz})', fontsize=13) - plt.tight_layout() - plt.savefig(os.path.join(out_dir,f'1d_{mode}_nx{nx}.png'), dpi=150); plt.close() - - # Coupes 2D - fig2, axes2 = plt.subplots(2, 3, figsize=(16,10)) - - sl_z = np.array([[nidx(i,j,mid_k) for i in range(nx)] for j in range(ny)]) - X2,Y2 = nodes[sl_z,0], nodes[sl_z,1] - for row,u,title in [(0,ux,'$u_x$ z=mid'),(1,uy,'$u_y$ z=mid')]: - im=axes2[row,0].pcolormesh(X2,Y2,u[sl_z],cmap='RdBu_r',shading='auto') - axes2[row,0].set(xlabel='x',ylabel='y',title=title) - plt.colorbar(im,ax=axes2[row,0]) - - sl_y = np.array([[nidx(i,mid_j,k) for i in range(nx)] for k in range(nz)]) - X2,Z2 = nodes[sl_y,0], nodes[sl_y,2] - for row,u,title in [(0,ux,'$u_x$ y=mid'),(1,uz,'$u_z$ y=mid')]: - im=axes2[row,1].pcolormesh(X2,Z2,u[sl_y],cmap='RdBu_r',shading='auto') - axes2[row,1].set(xlabel='x',ylabel='z',title=title) - plt.colorbar(im,ax=axes2[row,1]) - - sl_x = np.array([[nidx(mid_i,j,k) for j in range(ny)] for k in range(nz)]) - Y2,Z2 = nodes[sl_x,1], nodes[sl_x,2] - for row,u,title in [(0,uy,'$u_y$ x=mid'),(1,uz,'$u_z$ x=mid')]: - im=axes2[row,2].pcolormesh(Y2,Z2,u[sl_x],cmap='RdBu_r',shading='auto') - axes2[row,2].set(xlabel='y',ylabel='z',title=title) - plt.colorbar(im,ax=axes2[row,2]) - - plt.suptitle(f'MMS 3D Hexa Q1 [{mode}] — Displacement ({nx}×{ny}×{nz})', fontsize=13) - plt.tight_layout() - plt.savefig(os.path.join(out_dir,f'2d_{mode}_nx{nx}.png'), dpi=150); plt.close() - - - - - -def run_convergence_study(mode, nx_list=None, L=1.0, E=1e6, nu=0.3, out_dir=None): - - if nx_list is None: - nx_list = [3, 4, 6, 8, 11, 16] - - if out_dir is None: - out_dir = RESULTS_DIR - - os.makedirs(out_dir, exist_ok=True) - - hs, l2s, h1s = [], [], [] - - - for nx in nx_list: - ny = nz = nx - h = L / (nx - 1) - - print(f"[{mode}] Simulation nx={nx} (h={h:.4f})...") - - # Simulation - nodes, ux, uy, uz = run_simulation_3d(L, E, nu, nx, ny, nz, mode=mode) - - # ============ errors ============== - l2_err = l2_error_3d(nodes, ux, uy, uz, L, mode) - h1_err = h1_error_3d(nodes, ux, uy, uz, L, nx, ny, nz, mode) - - hs.append(h) - l2s.append(l2_err) - h1s.append(h1_err) - - if len(hs) == 1: - print(f" L2 error: {l2_err:.6e}") - print(f" H1 error: {h1_err:.6e}") - else: - l2_rate = np.log(l2s[-2]/l2s[-1]) / np.log(hs[-2]/hs[-1]) - h1_rate = np.log(h1s[-2]/h1s[-1]) / np.log(hs[-2]/hs[-1]) - print(f" L2 error: {l2_err:.6e} (rate {l2_rate:.2f})") - print(f" H1 error: {h1_err:.6e} (rate {h1_rate:.2f})") - - if nx == nx_list[-1]: - plot_displacement(nodes, ux, uy, uz, L, nx, ny, nz, mode, out_dir) - - plot_convergence(hs, l2s, h1s, mode, out_dir) - - - hs_arr = np.array(hs) - l2_arr = np.array(l2s) - h1_arr = np.array(h1s) - - l2_orders = np.diff(np.log(l2_arr)) / np.diff(np.log(hs_arr)) - h1_orders = np.diff(np.log(h1_arr)) / np.diff(np.log(hs_arr)) - - print(f"\n L2 : ordre moyen = {l2_orders.mean():.2f} (attendu: 2.0)") - print(f" H1 : ordre moyen = {h1_orders.mean():.2f} (attendu: 1.0)\n") - - return {'hs': hs, 'l2s': l2s, 'h1s': h1s} -