[MMS] Split 3D files in many (like 1D/2D) and reuse top level scripts

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]

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]
+        }
     }
 }

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/<stem>.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/<stem>.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)