[MMS] Unifying 1D to 3D refactor

examples/Freefem/MMS/1D/bar.py

@@ -20,6 +20,8 @@
     H1_QUADRATURE_1D,
 )
 from bar_solution import BarSolution1D
+from output import write_solution_table
+from scene import NodalForceAssembler
 
 RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results")
 
@@ -42,23 +44,13 @@ def load_params(path=None):
 # SOFA runner
 # ---------------------------------------------------------------------------
 
-class BodyForceAssembler(Sofa.Core.Controller):
-    """Fill the BodyForce field after init from nodes/edges read off the topology."""
-
-    def __init__(self, dofs, topology, body_force, f_body, quadrature, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-        self.dofs       = dofs
-        self.topology   = topology
-        self.body_force = body_force
-        self.f_body     = f_body
-        self.quadrature = quadrature
-
-    def onSimulationInitDoneEvent(self, event):
-        nodes = self.dofs.rest_position.array().copy().flatten()
-        edges = self.topology.edges.array().copy()
-        nodal_forces = assemble_nodal_forces_1d(self.f_body, nodes, edges, self.quadrature)
-        with self.body_force.forces.writeableArray() as forces:
-            forces[:, 0] = nodal_forces
+def _bar_force_compute(f_body, quadrature):
+    """Return a `(nodes, topology) -> (nx, 1) force array` for the assembler."""
+    def compute(nodes, topology):
+        flat  = nodes.copy().flatten()
+        edges = topology.edges.array().copy()
+        return assemble_nodal_forces_1d(f_body, flat, edges, quadrature).reshape(-1, 1)
+    return compute
 
 
 def build_bar_scene(root, mms, E_eff, nx):
@@ -122,12 +114,13 @@ def build_bar_scene(root, mms, E_eff, nx):
                                indices=list(range(nx)),
                                forces=[[0.0]] * nx)
 
-    Bar.addObject(BodyForceAssembler(
+    Bar.addObject(NodalForceAssembler(
         dofs=dofs,
         topology=Bar.topology,
-        body_force=body_force,
-        f_body=lambda xi: mms.source(xi, E_eff),
-        quadrature=mms.source_quadrature,
+        force_field=body_force,
+        compute_forces=_bar_force_compute(
+            f_body=lambda xi: mms.source(xi, E_eff),
+            quadrature=mms.source_quadrature),
         name="bodyForceAssembler"))
 
     mms.apply_bcs(Bar, E_eff, nx)
@@ -160,26 +153,6 @@ def solve_bar(mms, E_eff, nx):
 # Output helpers
 # ---------------------------------------------------------------------------
 
-def write_solution_table(case, x, u_h, u_ex, error_dict):
-    """
-    Write per-node solution table and error summary to results/<case>_solution.txt.
-
-    u_ex : callable x -> float
-    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"{case}_solution.txt")
-    with open(path, "w") as f:
-        f.write(f"{'x':>10} | {'u_h':>15} | {'u_ex':>15} | {'error':>15}\n")
-        f.write("-" * 60 + "\n")
-        for xi, ui in zip(x, u_h):
-            ue  = u_ex(xi)
-            f.write(f"{xi:10.4f} | {ui:15.6e} | {ue:15.6e} | {abs(ui - ue):15.6e}\n")
-        f.write("\n")
-        for label, val in error_dict.items():
-            f.write(f"{label:12s} = {val:.6e}\n")
-
-
 def plot_solution(case, x, u_h, u_ex, label_ex):
     """Save solution comparison plot to results/<case>_solution.png.
 
@@ -210,5 +183,6 @@ def run_reference_scene(mms):
     sol = solve_bar(mms, cfg["E_eff"], cfg["reference"]["nx"])
     l2  = l2_error_1d(sol.x0, sol.edges, sol.u_h, mms.u_ex, L2_QUADRATURE_1D)
     h1  = h1_semi_error_1d(sol.x0, sol.edges, sol.u_h, mms.du_ex, H1_QUADRATURE_1D)
-    write_solution_table(mms.name, sol.x0, sol.u_h, mms.u_ex, {"L2": l2, "H1_semi": h1})
+    write_solution_table(f"solution_{mms.name}", sol.x0, sol.u_h, mms.u_ex,
+                         RESULTS_DIR, {"L2": l2, "H1_semi": h1})
     plot_solution(mms.name, sol.x0, sol.u_h, mms.u_ex, mms.plot_label)

examples/Freefem/MMS/1D/manufactured_solution.py

@@ -1,11 +1,14 @@
 """Abstract base class for 1D MMS cases (non-dimensional, x ∈ [0,1])."""
 
-from abc import ABC, abstractmethod
+import os
+import sys
+from abc import abstractmethod
 
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from mms_case import MMSCase
 
-class MMSCase1D(ABC):
-    name              = None  # case identifier (must match the params.json key)
-    plot_label        = None  # LaTeX label for the exact solution
+
+class MMSCase1D(MMSCase):
     source_quadrature = None  # quadrature rule for body-force assembly only
 
     @abstractmethod

examples/Freefem/MMS/1D/run_convergence.py

@@ -1,13 +1,9 @@
 """Run the mesh-refinement convergence study for every 1D MMS case."""
 
-import os
-import numpy as np
-import matplotlib.pyplot as plt
-
-from quadratic  import mms as quadratic_mms
-from cubic      import mms as cubic_mms
-from sinusoidal import mms as sinusoidal_mms
-from exponential import mms as exponential_mms 
+from quadratic   import mms as quadratic_mms
+from cubic       import mms as cubic_mms
+from sinusoidal  import mms as sinusoidal_mms
+from exponential import mms as exponential_mms
 
 from bar import (
     RESULTS_DIR,
@@ -18,105 +14,29 @@
     L2_QUADRATURE_1D,
     H1_QUADRATURE_1D,
 )
-from plot_utils import annotate_convergence_rates
-
-
-def convergence_study(case, nx_list, run_fn, error_fns):
-    """
-    Run a mesh refinement convergence study and write results.
-
-    run_fn    : callable(nx) -> BarSolution1D
-    error_fns : dict { label: callable(x0, edges, u_h) -> float }
-    """
-    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── Convergence  {case} ──\n{hdr}", flush=True)
-
-    for k, nx_k in enumerate(nx_list):
-        h_k = 1.0 / (nx_k - 1)
-        sol = run_fn(nx_k)
-        row = {"nx": nx_k, "h": h_k}
-
-        for label, err_fn in error_fns.items():
-            e_k  = err_fn(sol.x0, sol.edges, sol.u_h)
-            rate = (f"{np.log(e_k / errors[label][-1]) / np.log(h_k / hs[-1]):.2f}"
-                    if k > 0 else "")
-            errors[label].append(e_k)
-            row[label]            = e_k
-            row[f"rate_{label}"]  = rate
-
-        line = f"{nx_k:5d} | {h_k: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_k)
-        rows.append(row)
-
-    stem = f"{case}_convergence"
-    write_convergence_table(stem, rows)
-    plot_series = [{"label": lbl, "errors": errors[lbl]} for lbl in err_labels]
-    plot_convergence(stem, hs, plot_series,
-                     title=f"Error Convergence — {case} function")
-
-
-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)
+from convergence import run_convergence_series
+from output      import plot_convergence
 
 
 if __name__ == "__main__":
     cfg = load_params()
     for mms in (quadratic_mms, cubic_mms, sinusoidal_mms, exponential_mms):
-        convergence_study(mms.name, cfg["convergence"]["nx_values"][mms.name],
-            run_fn    = lambda nx: solve_bar(mms, cfg["E_eff"], nx),
-            error_fns = {"L2":          lambda x, e, u: l2_error_1d(x, e, u, mms.u_ex, L2_QUADRATURE_1D),
-                         "H1 seminorm": lambda x, e, u: h1_semi_error_1d(x, e, u, mms.du_ex, H1_QUADRATURE_1D)})
+        stem = f"{mms.name}_convergence"
+        hs, errors = run_convergence_series(
+            nx_values  = cfg["convergence"]["nx_values"][mms.name],
+            run_fn     = lambda nx, _m=mms: solve_bar(_m, cfg["E_eff"], nx),
+            h_fn       = lambda nx: 1.0 / (nx - 1),
+            error_fns  = {
+                "L2": lambda sol, _m=mms: l2_error_1d(
+                    sol.x0, sol.edges, sol.u_h, _m.u_ex, L2_QUADRATURE_1D),
+                "H1": lambda sol, _m=mms: h1_semi_error_1d(
+                    sol.x0, sol.edges, sol.u_h, _m.du_ex, H1_QUADRATURE_1D),
+            },
+            banner     = f"-- Convergence  {mms.name} --",
+            results_dir = RESULTS_DIR,
+            table_stem  = stem,
+        )
+        plot_series = [{"label": lbl, "errors": errors[lbl]} for lbl in errors]
+        plot_convergence(stem, hs, plot_series,
+                         title=f"Error Convergence — {mms.name} function",
+                         results_dir=RESULTS_DIR)