[MMS] Refactor the 2D scripts

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

@@ -13,11 +13,11 @@
 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
 from fem import (
     line_quadrature,
-    assemble_nodal_forces,
-    l2_error,
-    h1_semi_error,
-    L2_QUADRATURE,
-    H1_QUADRATURE,
+    assemble_nodal_forces_1d,
+    l2_error_1d,
+    h1_semi_error_1d,
+    L2_QUADRATURE_1D,
+    H1_QUADRATURE_1D,
 )
 from bar_solution import BarSolution1D
 
@@ -56,7 +56,7 @@ def __init__(self, dofs, topology, body_force, f_body, quadrature, *args, **kwar
     def onSimulationInitDoneEvent(self, event):
         nodes = self.dofs.rest_position.array().copy().flatten()
         edges = self.topology.edges.array().copy()
-        nodal_forces = assemble_nodal_forces(self.f_body, nodes, edges, self.quadrature)
+        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
 
@@ -103,9 +103,9 @@ def build_bar_scene(root, mms, E_eff, nx):
                   newtonSolver="@newtonSolver",
                   linearSolver="@linearSolver")
 
-    Bar.addObject('EdgeSetTopologyContainer', name="topology"
-                  , edges="@../Grid/grid.edges"
-                  , position="@../Grid/grid.position")
+    Bar.addObject('EdgeSetTopologyContainer', name="topology",
+                  edges="@../Grid/grid.edges",
+                  position="@../Grid/grid.position")
 
     dofs = Bar.addObject('MechanicalObject',
                          name="dofs",
@@ -130,18 +130,14 @@ def build_bar_scene(root, mms, E_eff, nx):
         quadrature=mms.source_quadrature,
         name="bodyForceAssembler"))
 
-    Bar.addObject('FixedProjectiveConstraint', indices=0)
-    Bar.addObject('ConstantForceField',
-                  name="NeumannTip",
-                  indices=nx - 1,
-                  forces=mms.traction_bc(E_eff))
+    mms.apply_bcs(Bar, E_eff, nx)
 
 
 def case_scene(mms):
     """Return a `createScene(rootNode)` bound to this MMS case."""
     def createScene(rootNode):
         cfg = load_params()
-        build_bar_scene(rootNode, mms, cfg["E_eff"], cfg["nx"])
+        build_bar_scene(rootNode, mms, cfg["E_eff"], cfg["reference"]["nx"])
         return rootNode
     return createScene
 
@@ -211,8 +207,8 @@ def plot_solution(case, x, u_h, u_ex, label_ex):
 def run_reference_scene(mms):
     """Solve one MMS case at the reference mesh, write the solution table and plot."""
     cfg = load_params()
-    sol = solve_bar(mms, cfg["E_eff"], cfg["nx"])
-    l2  = l2_error(sol.x0, sol.edges, sol.u_h, mms.u_ex, L2_QUADRATURE)
-    h1  = h1_semi_error(sol.x0, sol.edges, sol.u_h, mms.du_ex, H1_QUADRATURE)
+    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})
     plot_solution(mms.name, sol.x0, sol.u_h, mms.u_ex, mms.plot_label)

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

@@ -27,6 +27,13 @@ def du_ex(self, xi):
     def source(self, xi, E):
         return E * (6.0 * xi - 2.0)
 
+    def apply_bcs(self, Bar, E_eff, nx):
+        Bar.addObject("FixedProjectiveConstraint", indices=0)
+        Bar.addObject("ConstantForceField",
+                      name="NeumannTip",
+                      indices=nx - 1,
+                      forces=self.traction_bc(E_eff))
+
 
 mms = Cubic()
 createScene = case_scene(mms)

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

@@ -26,6 +26,13 @@ def du_ex(self, xi):
     def source(self, xi, E):
         return -E * np.exp(xi)   # f = -E·u'' = -E·exp(x)
 
+    def apply_bcs(self, Bar, E_eff, nx):
+        Bar.addObject("FixedProjectiveConstraint", indices=0)
+        Bar.addObject("ConstantForceField",
+                      name="NeumannTip",
+                      indices=nx - 1,
+                      forces=self.traction_bc(E_eff))
+
 mms = Exponential()
 createScene = case_scene(mms)
 

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

@@ -23,3 +23,12 @@ def source(self, xi, E):
     def traction_bc(self, E):
         """Neumann traction at x=1: F_N = E · u'(1)."""
         return E * self.du_ex(1.0)
+
+    @abstractmethod
+    def apply_bcs(self, Bar, E_eff, nx):
+        """Install Dirichlet + Neumann BCs on the SOFA `Bar` node.
+
+        The MMS author chooses the BC pattern matching the manufactured
+        solution. Typical: Dirichlet at node 0, Neumann tip force at node
+        nx-1 using `self.traction_bc(E_eff)`.
+        """

examples/Freefem/MMS/1D/params.json

@@ -1,11 +1,15 @@
 {
     "length": 2.0,
     "youngModulus": 1e6,
-    "nx": 10,
-    "nxConvergence": {
-        "quadratic":   [2, 4, 6, 8, 10, 12, 16],
-        "cubic":       [2, 4, 6, 8, 10, 12, 16],
-        "sinusoidal":  [2, 4, 6, 8, 10, 12, 16],
-        "exponential" :[2, 4, 6, 8, 10, 12, 16]
+    "reference": {
+        "nx": 10
+    },
+    "convergence": {
+        "nx_values": {
+            "quadratic":   [2, 4, 6, 8, 10, 12, 16, 32],
+            "cubic":       [2, 4, 6, 8, 10, 12, 16, 32],
+            "sinusoidal":  [2, 4, 6, 8, 10, 12, 16, 32],
+            "exponential": [2, 4, 6, 8, 10, 12, 16, 32]
+        }
     }
 }

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

@@ -27,6 +27,13 @@ def du_ex(self, xi):
     def source(self, xi, E):
         return 2.0 * E
 
+    def apply_bcs(self, Bar, E_eff, nx):
+        Bar.addObject("FixedProjectiveConstraint", indices=0)
+        Bar.addObject("ConstantForceField",
+                      name="NeumannTip",
+                      indices=nx - 1,
+                      forces=self.traction_bc(E_eff))
+
 
 mms = Quadratic()
 createScene = case_scene(mms)

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

@@ -13,11 +13,12 @@
     RESULTS_DIR,
     load_params,
     solve_bar,
-    l2_error,
-    h1_semi_error,
-    L2_QUADRATURE,
-    H1_QUADRATURE,
+    l2_error_1d,
+    h1_semi_error_1d,
+    L2_QUADRATURE_1D,
+    H1_QUADRATURE_1D,
 )
+from plot_utils import annotate_convergence_rates
 
 
 def convergence_study(case, nx_list, run_fn, error_fns):
@@ -31,6 +32,11 @@ def convergence_study(case, nx_list, run_fn, error_fns):
     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)
@@ -44,22 +50,30 @@ def convergence_study(case, nx_list, run_fn, error_fns):
             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)
 
-    write_convergence_table(case, rows)
-    plot_convergence(case, hs, errors)
+    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(case, rows):
+def write_convergence_table(stem, rows):
     """
-    Write convergence table to results/<case>_convergence.txt.
+    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"{case}_convergence.txt")
+    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:
@@ -73,41 +87,36 @@ def write_convergence_table(case, rows):
             f.write(line + "\n")
 
 
-def plot_convergence(case, hs, error_series):
+def plot_convergence(stem, hs, series, title, ylabel="Error"):
     """
-    Save log-log convergence plot to results/<case>_convergence.png.
+    Save log-log convergence plot to results/<stem>.png.
 
-    error_series : dict { label: array-like of errors }
+    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)
-    markers = ["bo-", "rs--", "g^:"]
+    default = ["bo-", "rs--", "g^:", "m^-"]
     fig, ax = plt.subplots(figsize=(8, 5))
-    for (label, errors), marker in zip(error_series.items(), markers):
-        e_arr = np.array(errors)
-        ax.loglog(h_arr, e_arr, marker, label=label, linewidth=2, markersize=7)
-        for k in range(1, len(h_arr)):
-            rate  = np.log(e_arr[k] / e_arr[k-1]) / np.log(h_arr[k] / h_arr[k-1])
-            x_mid = np.sqrt(h_arr[k] * h_arr[k-1])
-            y_mid = np.sqrt(e_arr[k] * e_arr[k-1])
-            ax.annotate(f"{rate:.2f}", xy=(x_mid, y_mid),
-                        xytext=(0, 8), textcoords="offset points",
-                        ha="center", fontsize=9)
+    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("Error")
-    ax.set_title(f"Error Convergence — {case} function")
+    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"{case}_convergence.png"), dpi=150)
+    fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150)
     plt.close(fig)
 
 
 if __name__ == "__main__":
     cfg = load_params()
     for mms in (quadratic_mms, cubic_mms, sinusoidal_mms, exponential_mms):
-        convergence_study(mms.name, cfg["nxConvergence"][mms.name],
+        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(x, e, u, mms.u_ex, L2_QUADRATURE),
-                         "H1 seminorm": lambda x, e, u: h1_semi_error(x, e, u, mms.du_ex, H1_QUADRATURE)})
+            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)})

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

@@ -29,6 +29,13 @@ def du_ex(self, xi):
     def source(self, xi, E):
         return 4.0 * np.pi**2 * E * np.sin(2.0 * np.pi * xi)
 
+    def apply_bcs(self, Bar, E_eff, nx):
+        Bar.addObject("FixedProjectiveConstraint", indices=0)
+        Bar.addObject("ConstantForceField",
+                      name="NeumannTip",
+                      indices=nx - 1,
+                      forces=self.traction_bc(E_eff))
+
 
 mms = Sinusoidal()
 createScene = case_scene(mms)