diff --git a/examples/Freefem/MMS/1D/bar.py b/examples/Freefem/MMS/1D/bar.py index b1fe7b50..3904cb11 100644 --- a/examples/Freefem/MMS/1D/bar.py +++ b/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/_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/_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) diff --git a/examples/Freefem/MMS/1D/manufactured_solution.py b/examples/Freefem/MMS/1D/manufactured_solution.py index 1f8c9789..5366aeba 100644 --- a/examples/Freefem/MMS/1D/manufactured_solution.py +++ b/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 diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py index 0e3a51f0..d063d38d 100644 --- a/examples/Freefem/MMS/1D/run_convergence.py +++ b/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/.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) +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) diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py index 641f593b..420ccee4 100644 --- a/examples/Freefem/MMS/2D/beam.py +++ b/examples/Freefem/MMS/2D/beam.py @@ -11,15 +11,11 @@ # 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_2d, - assemble_traction_2d, - l2_error_2d, - h1_semi_error_2d, - quad_q1_rule, - tri_p1_rule, -) +from fem import quad_q1_rule, tri_p1_rule # re-exported for case files +from elements import element_quad, element_tri # re-exported for case files from beam_solution import BeamSolution2D +from output import write_solution_table +from scene import NodalForceAssembler RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results") @@ -52,173 +48,23 @@ def _dim_template(dim): return "Vec3d" if dim == "3d" else "Vec2d" -def _boundary_edges(nx, ny): - """Return (bottom, top, left, right) edge lists for a structured nx×ny grid.""" - bottom = [(i, i + 1) for i in range(nx - 1)] - top = [((ny - 1) * nx + i, (ny - 1) * nx + i + 1) for i in range(nx - 1)] - left = [(j * nx, (j + 1) * nx) for j in range(ny - 1)] - right = [(j * nx + (nx - 1), (j + 1) * nx + (nx - 1)) for j in range(ny - 1)] - return bottom, top, left, right - - -# --------------------------------------------------------------------------- -# 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. quad_q1_rule(2)) -# add_topology(rootNode, Beam, nx, ny, L) : 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. -# to_triangles(conn) : triangulation for matplotlib tricontourf. -# -# 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_2d, conn, mms, L, E, nu, nx, ny, dim): - xy = nodes_2d[:, :2] - - F = assemble_nodal_forces_2d( - lambda x, y: mms.source(x, y, E, nu, L, dim), - xy, conn, cls._source_rule(mms)) - - bottom, top, left, right = _boundary_edges(nx, ny) - sides = [(bottom, 0.0, -1.0), - (top, 0.0, +1.0), - (left, -1.0, 0.0), - (right, +1.0, 0.0)] - for edges, nrm_x, nrm_y in sides: - F += assemble_traction_2d( - lambda x, y, nx=nrm_x, ny=nrm_y: - mms.traction(x, y, nx, ny, E, nu, L, dim), - xy, edges) - return F - - @classmethod - def compute_l2(cls, sol, mms, L): - return l2_error_2d( - sol.nodes, sol.conn, np.column_stack([sol.ux, sol.uy]), - lambda x, y: mms.u_ex(x, y, L), - cls.ELEMENT_RULE) - - @classmethod - def compute_h1(cls, sol, mms, L): - return h1_semi_error_2d( - sol.nodes, sol.conn, np.column_stack([sol.ux, sol.uy]), - lambda x, y: mms.grad_u_ex(x, y, L), - cls.ELEMENT_RULE) - - -class _QuadElement(_ElementBase): - LABEL = "Q1 quad" - ELEMENT_RULE = staticmethod(quad_q1_rule(2)) # used for L²/H¹ error norms - - @staticmethod - def _source_rule(mms): - rule = mms.source_quadrature_quad - if rule is None: - raise ValueError( - f"{type(mms).__name__}.source_quadrature_quad must be set") - return rule - - @staticmethod - def add_topology(Beam): - topology = Beam.addObject("QuadSetTopologyContainer", name="topology", - quads="@../Grid/grid.quads", - position="@../Grid/grid.position") - Beam.addObject("QuadSetTopologyModifier") - return topology - - @staticmethod - def read_connectivity(topology): - return topology.quads.array().copy() - - @staticmethod - def to_triangles(conn): - tris = [] - for q in conn: - tris.append([q[0], q[1], q[2]]) - tris.append([q[0], q[2], q[3]]) - return np.array(tris) - - -class _TriElement(_ElementBase): - LABEL = "P1 tri" - ELEMENT_RULE = staticmethod(tri_p1_rule(3)) # used for L²/H¹ error norms - - @staticmethod - def _source_rule(mms): - rule = mms.source_quadrature_tri - if rule is None: - raise ValueError( - f"{type(mms).__name__}.source_quadrature_tri must be set") - return rule - - @staticmethod - def add_topology(Beam): - topology = Beam.addObject("TriangleSetTopologyContainer", name="topology") - Beam.addObject("Quad2TriangleTopologicalMapping", - input="@../Grid/grid", output="@topology") - Beam.addObject("TriangleSetTopologyModifier") - return topology - - @staticmethod - def read_connectivity(topology): - return topology.triangles.array().copy() - - @staticmethod - def to_triangles(conn): - return conn - - # --------------------------------------------------------------------------- -# Force assembly controller -# -# SOFA topology components (RegularGridTopology, QuadSetTopologyContainer, -# Quad2TriangleTopologicalMapping, ...) are only populated after init runs, -# not during Python scene-build time. This controller defers nodal-force -# assembly to onSimulationInitDoneEvent, where it reads positions off the -# MechanicalObject and connectivity off the SOFA topology, then fills the -# ConstantForceField that was added with placeholder zeros. -# -# Mirrors the 1D pattern in `bar.py:BodyForceAssembler`. +# SOFA scene # --------------------------------------------------------------------------- -class NodalForceAssembler(Sofa.Core.Controller): - def __init__(self, dofs, topology, force_field, element, mms, - L, E, nu, nx, ny, dim, *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.dim = nx, ny, dim - - def onSimulationInitDoneEvent(self, event): - nodes = self.dofs.rest_position.array().copy() - conn = self.element.read_connectivity(self.topology) - F_xy = self.element.compute_nodal_forces( - nodes, conn, self.mms, - self.L, self.E, self.nu, self.nx, self.ny, self.dim) - if self.dim == "3d": - F_full = np.hstack([F_xy, np.zeros((len(F_xy), 1))]) - else: - F_full = F_xy - with self.force_field.forces.writeableArray() as forces: - forces[:] = F_full - +def _beam_force_compute(element, mms, L, E, nu, nx, ny, dim): + """Return a `(nodes, topology) -> force array` for the shared NodalForceAssembler. -# --------------------------------------------------------------------------- -# SOFA scene -# --------------------------------------------------------------------------- + Pads the 2D force array with a zero z-column when the scene uses Vec3d. + """ + def compute(nodes, topology): + conn = element.read_connectivity(topology) + F_xy = element.compute_nodal_forces( + nodes, conn, mms, L, E, nu, nx, ny, dim) + if dim == "3d": + return np.hstack([F_xy, np.zeros((len(F_xy), 1))]) + return F_xy + return compute def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, nx=10, ny=10, with_visual=True, dim="2d"): @@ -294,8 +140,7 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, Beam.addObject(NodalForceAssembler( dofs=dofs, topology=topology, force_field=force_field, - element=element, mms=mms, - L=L, E=E, nu=nu, nx=nx, ny=ny, dim=dim, + compute_forces=_beam_force_compute(element, mms, L, E, nu, nx, ny, dim), name="nodalForceAssembler")) return dofs, topology @@ -329,29 +174,6 @@ def solve_beam(elem, mms, L, E, nu, nx, ny, dim="2d"): # Output helpers (mirror 1D bar.py) # --------------------------------------------------------------------------- -def write_solution_table(stem, x, y, ux_h, uy_h, u_ex, error_dict): - """ - Write per-node solution table and error summary to results/.txt. - - u_ex : callable (x, y) -> (ux_ex, uy_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} | {'ux_h':>15} | {'uy_h':>15} | " - f"{'ux_ex':>15} | {'uy_ex':>15} | {'err_x':>15} | {'err_y':>15}\n") - f.write("-" * 124 + "\n") - for xi, yi, uxi, uyi in zip(x, y, ux_h, uy_h): - uxe, uye = u_ex(xi, yi) - f.write(f"{xi:10.4f} | {yi:10.4f} | {uxi:15.6e} | {uyi:15.6e} | " - f"{uxe:15.6e} | {uye:15.6e} | " - f"{abs(uxi - uxe):15.6e} | {abs(uyi - uye):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, label, dim, hyp, nu, l2, h1): """Save 1-D mid-height profile (ux, uy vs x) to results/.png.""" os.makedirs(RESULTS_DIR, exist_ok=True) @@ -435,10 +257,10 @@ def run_reference_scene(elem, mms): stem = f"{mms.name}_{tag}_{dim}_nu{nu}_nx{nx}" xy = sol.nodes[:, :2] - write_solution_table(f"solution_{stem}", xy[:, 0], xy[:, 1], - sol.ux, sol.uy, + write_solution_table(f"solution_{stem}", xy, + np.column_stack([sol.ux, sol.uy]), lambda xi, yi: mms.u_ex(xi, yi, L), - {"L2": l2, "H1_semi": h1}) + RESULTS_DIR, {"L2": l2, "H1_semi": h1}) plot_solution_profile(f"solution_{stem}", sol, mms, L, nx, ny, label, dim, hyp, nu, l2, h1) plot_solution_fields (f"fields2D_{stem}", sol, elem, mms, L, nx, @@ -462,11 +284,3 @@ def createScene(rootNode): with_visual=True, dim=ref["dim"]) return rootNode return createScene - - -# ───────────────────────────────────────────────────────────────────────────── -# Element instances -# ───────────────────────────────────────────────────────────────────────────── - -element_quad = _QuadElement() -element_tri = _TriElement() diff --git a/examples/Freefem/MMS/2D/manufactured_solution.py b/examples/Freefem/MMS/2D/manufactured_solution.py index a03db6fe..2b90ca7e 100644 --- a/examples/Freefem/MMS/2D/manufactured_solution.py +++ b/examples/Freefem/MMS/2D/manufactured_solution.py @@ -1,6 +1,11 @@ """Abstract base class for 2D MMS cases (Cartesian domain [0,L]^2).""" -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 def lame(E, nu, dim): @@ -18,10 +23,7 @@ def lame(E, nu, dim): return lam, mu -class MMSCase2D(ABC): - name = None # case identifier (must match the params.json key) - plot_label = None # LaTeX label for the exact solution - +class MMSCase2D(MMSCase): # Body-force quadrature rules — must be set by each concrete case. # No framework fallback; assembly raises if either is unset. source_quadrature_quad = None # element rule for Q1 quads (e.g. quad_q1_rule(2)) diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py index f877c2c1..bad05b1f 100644 --- a/examples/Freefem/MMS/2D/run_convergence.py +++ b/examples/Freefem/MMS/2D/run_convergence.py @@ -4,10 +4,6 @@ text tables and convergence plots into the shared `results/` directory. """ -import os -import numpy as np -import matplotlib.pyplot as plt - from cubic import mms as cubic_mms from trigonometric import mms as trig_mms from incompressible import mms as incomp_mms @@ -19,12 +15,13 @@ element_quad, element_tri, ) -from plot_utils import annotate_convergence_rates +from convergence import run_convergence_series +from output import plot_convergence def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"): """ - Run convergence study for each element type in elem_specs, write a + Run a convergence series for each element type in elem_specs, write a per-(element) text table, and one shared plot with L²/H¹ for every element on the same axes. @@ -33,8 +30,6 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"): """ hyp = "plane strain" if dim == "3d" else "plane stress" print(f"\n PoissonRatio = {nu} ({hyp})", 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: @@ -42,86 +37,29 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"): tag = label.replace(" ", "_") stem = f"convergence_{mms.name}_{tag}_{dim}_nu{nu}" - print(f"\n── {label} [{dim} / {hyp}] {mms.name} nu={nu} ──\n{hdr}", - flush=True) - - rows, hs, l2s, h1s = [], [], [], [] - for k, nx in enumerate(nx_values): - ny = nx - h = L / (nx - 1) - sol = solve_beam(elem, mms, L, E, nu, nx, ny, dim=dim) - l2 = elem.compute_l2(sol, mms, L) - h1 = elem.compute_h1(sol, mms, L) + hs, errors = run_convergence_series( + nx_values = nx_values, + run_fn = lambda nx, _e=elem: solve_beam( + _e, mms, L, E, nu, nx, nx, dim=dim), + h_fn = lambda nx: L / (nx - 1), + error_fns = { + "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L), + "H1": lambda sol, _e=elem: _e.compute_h1(sol, mms, L), + }, + banner = f"-- {label} [{dim} / {hyp}] {mms.name} nu={nu} --", + results_dir = RESULTS_DIR, + table_stem = stem, + ) - 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"]}) + "errors": errors["L2"], "style": spec["l2_style"]}) plot_series.append({"label": f"{label} H¹", - "errors": h1s, "style": spec["h1_style"]}) + "errors": errors["H1"], "style": spec["h1_style"]}) hs_ref = hs title = f"Convergence — {mms.name} [{dim} / {hyp}] nu={nu}" plot_convergence(f"convergence_{mms.name}_{dim}_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) + hs_ref, plot_series, title=title, results_dir=RESULTS_DIR) if __name__ == "__main__": @@ -140,7 +78,7 @@ def plot_convergence(stem, hs, series, title, ylabel="Error"): # dim="2d" → Vec2d / plane stress; dim="3d" → Vec3d / plane strain for mms in (cubic_mms, trig_mms, incomp_mms): nx_vals = conv["nx_values"][mms.name] - print(f"\n══ {mms.name} ══") + print(f"\n== {mms.name} ==") for DIM in conv["dim_values"]: for nu in conv["nu_values"]: convergence_study(specs, mms, L, E, nu, nx_vals, dim=DIM) diff --git a/examples/Freefem/MMS/3D/manufactured_solution.py b/examples/Freefem/MMS/3D/manufactured_solution.py index de39c6d7..6d64839f 100644 --- a/examples/Freefem/MMS/3D/manufactured_solution.py +++ b/examples/Freefem/MMS/3D/manufactured_solution.py @@ -1,8 +1,13 @@ """Abstract base class for 3D MMS cases (Cartesian domain [0,L]^3).""" -from abc import ABC, abstractmethod +import os +import sys +from abc import abstractmethod import numpy as np +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +from mms_case import MMSCase + def lame(E, nu): """Lamé parameters (lambda, mu) for 3D linear elasticity (full Hooke).""" @@ -11,10 +16,7 @@ def lame(E, 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 - +class MMSCase3D(MMSCase): # Body-force quadrature rule — must be set by each concrete case. # No framework fallback; assembly raises if unset. source_quadrature_hex = None # element rule for Q1 hexes (e.g. hex_q1_rule(2)) diff --git a/examples/Freefem/MMS/3D/run_convergence.py b/examples/Freefem/MMS/3D/run_convergence.py index 84b9dee5..44591a83 100644 --- a/examples/Freefem/MMS/3D/run_convergence.py +++ b/examples/Freefem/MMS/3D/run_convergence.py @@ -5,10 +5,6 @@ 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 ( @@ -17,20 +13,19 @@ solve_solid, element_hex, ) -from plot_utils import annotate_convergence_rates +from convergence import run_convergence_series +from output import plot_convergence def convergence_study(elem_specs, mms, L, E, nu, nx_values): """ - Run convergence study for each element type in elem_specs, write a + Run a convergence series for each element type in elem_specs, write a per-(element) text table, and one shared plot with L²/H¹ for every element on the same axes. elem_specs : list of dicts with keys 'elem', 'label', 'l2_style', 'h1_style' """ 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: @@ -38,85 +33,29 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values): 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) + hs, errors = run_convergence_series( + nx_values = nx_values, + run_fn = lambda nx, _e=elem: solve_solid( + _e, mms, L, E, nu, nx, nx, nx), + h_fn = lambda nx: L / (nx - 1), + error_fns = { + "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L), + "H1": lambda sol, _e=elem: _e.compute_h1(sol, mms, L), + }, + banner = f"-- {label} {mms.name} nu={nu} --", + results_dir = RESULTS_DIR, + table_stem = stem, + ) - 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"]}) + "errors": errors["L2"], "style": spec["l2_style"]}) plot_series.append({"label": f"{label} H¹", - "errors": h1s, "style": spec["h1_style"]}) + "errors": errors["H1"], "style": spec["h1_style"]}) hs_ref = hs title = f"Convergence — {mms.name} nu={nu}" plot_convergence(f"convergence_{mms.name}_nu{nu}", - hs_ref, plot_series, title=title) - - -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) + hs_ref, plot_series, title=title, results_dir=RESULTS_DIR) if __name__ == "__main__": @@ -132,6 +71,6 @@ def plot_convergence(stem, hs, series, title, ylabel="Error"): for mms in (sinus_neumann_mms,): nx_vals = conv["nx_values"][mms.name] - print(f"\n══ {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/solid.py b/examples/Freefem/MMS/3D/solid.py index ce403450..7fb68bec 100644 --- a/examples/Freefem/MMS/3D/solid.py +++ b/examples/Freefem/MMS/3D/solid.py @@ -11,14 +11,11 @@ # 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 fem import hex_q1_rule # re-exported for case files +from elements import element_hex # re-exported for case files from solid_solution import SolidSolution3D +from output import write_solution_table +from scene import NodalForceAssembler RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results") @@ -46,151 +43,18 @@ def get_nodes_3d(L, nx, ny, nz): 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 _solid_force_compute(element, mms, L, E, nu, nx, ny, nz): + """Return a `(nodes, topology) -> (N, 3) force array` for NodalForceAssembler.""" + def compute(nodes, topology): + conn = element.read_connectivity(topology) + return element.compute_nodal_forces( + nodes, conn, mms, L, E, nu, nx, ny, nz) + return compute + def build_solid_scene(rootNode, mms, element, 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. @@ -256,8 +120,7 @@ def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, 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, + compute_forces=_solid_force_compute(element, mms, L, E, nu, nx, ny, nz), name="nodalForceAssembler")) return dofs, topology @@ -291,33 +154,6 @@ def solve_solid(elem, mms, L, E, nu, nx, ny, nz): # 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) @@ -439,11 +275,10 @@ def run_reference_scene(elem, mms): 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, + write_solution_table(f"solution_{stem}", xyz, + np.column_stack([sol.ux, sol.uy, sol.uz]), lambda xi, yi, zi: mms.u_ex(xi, yi, zi, L), - {"L2": l2, "H1_semi": h1}) + RESULTS_DIR, {"L2": l2, "H1_semi": h1}) plot_solution_profile(f"solution_{stem}", sol, mms, L, nx, ny, nz, label, nu, l2, h1) plot_solution_slices (f"fields3D_{stem}", sol, mms, L, nx, ny, nz, @@ -468,10 +303,3 @@ def createScene(rootNode): with_visual=True) return rootNode return createScene - - -# ───────────────────────────────────────────────────────────────────────────── -# Element instances -# ───────────────────────────────────────────────────────────────────────────── - -element_hex = _HexElement() diff --git a/examples/Freefem/MMS/convergence.py b/examples/Freefem/MMS/convergence.py new file mode 100644 index 00000000..ddf98824 --- /dev/null +++ b/examples/Freefem/MMS/convergence.py @@ -0,0 +1,62 @@ +"""Shared convergence-study helpers for MMS 1D/2D/3D drivers.""" + +import numpy as np + +from output import write_convergence_table + + +def run_convergence_series(nx_values, run_fn, h_fn, error_fns, banner, + results_dir, table_stem): + """ + Run a mesh-refinement series and write a table. + + For each `nx` in `nx_values`, runs `run_fn(nx)` to get a solution, + evaluates every error in `error_fns`, prints a live row, and writes + the assembled table to ``/.txt``. + + Parameters + ---------- + nx_values : iterable of mesh-refinement parameters + run_fn : callable(nx) -> Solution (anything `error_fns` understands) + h_fn : callable(nx) -> mesh size h (drives the rate denominator + and the printed h column) + error_fns : ordered dict { label: callable(sol) -> float } + banner : string printed before the table header + results_dir, table_stem : the file ``/.txt`` + receives the assembled table + + Returns + ------- + (hs, errors) : list of mesh sizes; dict {label: list of errors}. + """ + hs = [] + errors = {label: [] for label in error_fns} + rows = [] + err_labels = list(error_fns.keys()) + + hdr = f"{'nx':>5} | {'h':>10}" + "".join( + f" | {lbl:>14} | {('rate_' + lbl):>7}" for lbl in err_labels) + print(f"\n{banner}\n{hdr}", flush=True) + + for k, nx in enumerate(nx_values): + h = h_fn(nx) + sol = run_fn(nx) + row = {"nx": nx, "h": h} + for label, err_fn in error_fns.items(): + e_k = err_fn(sol) + rate = (f"{np.log(e_k / errors[label][-1]) / np.log(h / hs[-1]):.2f}" + if k > 0 else "") + errors[label].append(e_k) + row[label] = e_k + row[f"rate_{label}"] = rate + + line = f"{nx:5d} | {h:10.4f}" + for lbl in err_labels: + line += f" | {row[lbl]:14.6e} | {row[f'rate_{lbl}']:>7}" + print(line, flush=True) + + hs.append(h) + rows.append(row) + + write_convergence_table(table_stem, rows, results_dir) + return hs, errors diff --git a/examples/Freefem/MMS/elements.py b/examples/Freefem/MMS/elements.py new file mode 100644 index 00000000..1f4fc9c4 --- /dev/null +++ b/examples/Freefem/MMS/elements.py @@ -0,0 +1,259 @@ +"""Element strategies for MMS scenes. + +Each strategy declares: + + LABEL : human-readable identifier (used in plot legends) + ELEMENT_RULE : an element rule from `fem.py` driving the L²/H¹ norms + _source_rule(mms) : the body-force quadrature rule (read off the mms; + raises if the case did not set it) + add_topology(parent_node) : add the SOFA topology container and return it + read_connectivity(topology) : extract the connectivity array post-init + compute_nodal_forces(...) : assemble body + traction nodal forces + compute_l2(sol, mms, L) : L² error norm on a solution + compute_h1(sol, mms, L) : H¹ semi-norm error on a solution + +Two `_ElementBase` classes live here — one for 2D (4 boundary edges, +dim-dependent plane stress / plane strain) and one for 3D (6 boundary +quads, single constitutive branch). Concrete elements (`_QuadElement`, +`_TriElement`, `_HexElement`) subclass the matching base and pin the +topology container choice and the rules. +""" + +import numpy as np + +from fem import ( + assemble_nodal_forces, + assemble_traction, + l2_error, + h1_semi_error, + quad_q1_rule, + tri_p1_rule, + hex_q1_rule, + edge_line_rule, + quad_face_rule, +) + + +# --------------------------------------------------------------------------- +# Boundary-facet helpers +# --------------------------------------------------------------------------- + +def _boundary_edges(nx, ny): + """Return (bottom, top, left, right) edge lists for a structured nx×ny grid.""" + bottom = [(i, i + 1) for i in range(nx - 1)] + top = [((ny - 1) * nx + i, (ny - 1) * nx + i + 1) for i in range(nx - 1)] + left = [(j * nx, (j + 1) * nx) for j in range(ny - 1)] + right = [(j * nx + (nx - 1), (j + 1) * nx + (nx - 1)) for j in range(ny - 1)] + return bottom, top, left, right + + +def _boundary_quads(nx, ny, nz): + """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 using SOFA's regular-grid index + convention `idx(i, j, k) = i + j*nx + k*nx*ny`. Orientation is consistent + per face but its sign does not affect `assemble_traction`, 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 + + +# --------------------------------------------------------------------------- +# 2D elements (Vec2d / plane stress OR Vec3d / plane strain via dim arg) +# --------------------------------------------------------------------------- + +class _ElementBase2D: + @classmethod + def compute_nodal_forces(cls, nodes_2d, conn, mms, L, E, nu, nx, ny, dim): + xy = nodes_2d[:, :2] + + F = assemble_nodal_forces( + lambda x, y: mms.source(x, y, E, nu, L, dim), + xy, conn, cls._source_rule(mms)) + + bottom, top, left, right = _boundary_edges(nx, ny) + sides = [(bottom, 0.0, -1.0), + (top, 0.0, +1.0), + (left, -1.0, 0.0), + (right, +1.0, 0.0)] + edge_rule = edge_line_rule(2) + for edges, nrm_x, nrm_y in sides: + F += assemble_traction( + lambda x, y, nx=nrm_x, ny=nrm_y: + mms.traction(x, y, nx, ny, E, nu, L, dim), + xy, edges, edge_rule) + return F + + @classmethod + def compute_l2(cls, sol, mms, L): + xy = sol.nodes[:, :2] + return l2_error( + xy, sol.conn, np.column_stack([sol.ux, sol.uy]), + lambda x, y: mms.u_ex(x, y, L), + cls.ELEMENT_RULE) + + @classmethod + def compute_h1(cls, sol, mms, L): + xy = sol.nodes[:, :2] + return h1_semi_error( + xy, sol.conn, np.column_stack([sol.ux, sol.uy]), + lambda x, y: mms.grad_u_ex(x, y, L), + cls.ELEMENT_RULE) + + +class _QuadElement(_ElementBase2D): + LABEL = "Q1 quad" + ELEMENT_RULE = staticmethod(quad_q1_rule(2)) # used for L²/H¹ error norms + + @staticmethod + def _source_rule(mms): + rule = mms.source_quadrature_quad + if rule is None: + raise ValueError( + f"{type(mms).__name__}.source_quadrature_quad must be set") + return rule + + @staticmethod + def add_topology(Beam): + topology = Beam.addObject("QuadSetTopologyContainer", name="topology", + quads="@../Grid/grid.quads", + position="@../Grid/grid.position") + Beam.addObject("QuadSetTopologyModifier") + return topology + + @staticmethod + def read_connectivity(topology): + return topology.quads.array().copy() + + @staticmethod + def to_triangles(conn): + tris = [] + for q in conn: + tris.append([q[0], q[1], q[2]]) + tris.append([q[0], q[2], q[3]]) + return np.array(tris) + + +class _TriElement(_ElementBase2D): + LABEL = "P1 tri" + ELEMENT_RULE = staticmethod(tri_p1_rule(3)) # used for L²/H¹ error norms + + @staticmethod + def _source_rule(mms): + rule = mms.source_quadrature_tri + if rule is None: + raise ValueError( + f"{type(mms).__name__}.source_quadrature_tri must be set") + return rule + + @staticmethod + def add_topology(Beam): + topology = Beam.addObject("TriangleSetTopologyContainer", name="topology") + Beam.addObject("Quad2TriangleTopologicalMapping", + input="@../Grid/grid", output="@topology") + Beam.addObject("TriangleSetTopologyModifier") + return topology + + @staticmethod + def read_connectivity(topology): + return topology.triangles.array().copy() + + @staticmethod + def to_triangles(conn): + return conn + + +# --------------------------------------------------------------------------- +# 3D elements (Vec3d / full Hooke) +# --------------------------------------------------------------------------- + +class _ElementBase3D: + @classmethod + def compute_nodal_forces(cls, nodes_3d, conn, mms, L, E, nu, nx, ny, nz): + xyz = nodes_3d[:, :3] + + F = assemble_nodal_forces( + lambda x, y, z: mms.source(x, y, z, E, nu, L), + xyz, conn, cls._source_rule(mms)) + + xm, xp, ym, yp, zm, zp = _boundary_quads(nx, ny, nz) + sides = [(xm, -1.0, 0.0, 0.0), + (xp, +1.0, 0.0, 0.0), + (ym, 0.0, -1.0, 0.0), + (yp, 0.0, +1.0, 0.0), + (zm, 0.0, 0.0, -1.0), + (zp, 0.0, 0.0, +1.0)] + face_rule = quad_face_rule(2) + for quads, nrm_x, nrm_y, nrm_z in sides: + F += assemble_traction( + lambda x, y, z, nx=nrm_x, ny=nrm_y, nz=nrm_z: + mms.traction(x, y, z, nx, ny, nz, E, nu, L), + xyz, quads, face_rule) + return F + + @classmethod + def compute_l2(cls, sol, mms, L): + return l2_error( + sol.nodes, sol.conn, + np.column_stack([sol.ux, sol.uy, sol.uz]), + lambda x, y, z: mms.u_ex(x, y, z, L), + cls.ELEMENT_RULE) + + @classmethod + def compute_h1(cls, sol, mms, L): + return h1_semi_error( + sol.nodes, sol.conn, + np.column_stack([sol.ux, sol.uy, sol.uz]), + lambda x, y, z: mms.grad_u_ex(x, y, z, L), + cls.ELEMENT_RULE) + + +class _HexElement(_ElementBase3D): + LABEL = "Q1 hex" + ELEMENT_RULE = staticmethod(hex_q1_rule(2)) # 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() + + +# --------------------------------------------------------------------------- +# Instances (one per element type) +# --------------------------------------------------------------------------- + +element_quad = _QuadElement() +element_tri = _TriElement() +element_hex = _HexElement() diff --git a/examples/Freefem/MMS/fem.py b/examples/Freefem/MMS/fem.py index 8288d724..edb7e0b1 100644 --- a/examples/Freefem/MMS/fem.py +++ b/examples/Freefem/MMS/fem.py @@ -1,12 +1,38 @@ """FEM toolbox: quadrature, assembly, error norms. + +Two parallel families of routines: + +* **1D scalar-integrand** path used by the 1D bar driver: + `line_quadrature`, `assemble_nodal_forces_1d`, `l2_error_1d`, + `h1_semi_error_1d`. + +* **Dim-agnostic vector path** used by the 2D and 3D drivers: + `quad_q1_rule`, `tri_p1_rule`, `hex_q1_rule` (element rules), + `edge_line_rule`, `quad_face_rule` (facet rules), + `assemble_nodal_forces`, `assemble_traction`, + `l2_error`, `h1_semi_error`. + +Element-rule protocol: + + rule(xe) yields (coords, w, N, dN_phys) + xe : (n_local, dim) physical node coordinates of one element + coords : (dim,) physical Gauss-point coordinates + w : scalar weight (already × detJ / triangle area) + N : (n_local,) shape values at the Gauss point + dN_phys : (dim, n_local) physical-coord shape gradients, indexed + by (spatial axis, local node) + +Facet-rule protocol: + + rule(xe) yields (coords, w, N) + same convention; no shape gradient (not needed on boundary). """ import numpy as np # --------------------------------------------------------------------------- -# Quadrature rules +# Canonical 1D Gauss-Legendre table on [-1, 1] # --------------------------------------------------------------------------- -# Canonical Gauss-Legendre points and weights on the reference segment [-1, 1]. _GAUSS_LEGENDRE_1D = { 1: (np.array([0.0]), np.array([2.0])), @@ -15,6 +41,10 @@ } +# --------------------------------------------------------------------------- +# 1D scalar-integrand quadrature (consumed by the 1D bar driver only) +# --------------------------------------------------------------------------- + def line_quadrature(n_pts): """Return a Gauss-Legendre rule on a physical segment using `n_pts` points. @@ -38,11 +68,11 @@ def rule(g, x1, x2): # --------------------------------------------------------------------------- -# 2D Q1 reference shape functions on [-1,1]^2. -# Node ordering: (-1,-1), (+1,-1), (+1,+1), (-1,+1) +# Reference-element shape functions # --------------------------------------------------------------------------- def _shape_q1(xi, eta): + """2D Q1 shape on [-1, 1]^2. Node order: (-,-), (+,-), (+,+), (-,+).""" N = 0.25 * np.array([ (1 - xi) * (1 - eta), (1 + xi) * (1 - eta), @@ -50,92 +80,14 @@ def _shape_q1(xi, 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)]) + dN_deta = 0.25 * np.array([-(1 - xi), -(1 + xi), (1 + xi), (1 - xi) ]) return N, dN_dxi, dN_deta -# --------------------------------------------------------------------------- -# 2D element rules -# -# An "element rule" is a callable rule(xe, ye) that yields one tuple -# (xg, yg, w, N, dN_dx, dN_dy) -# per Gauss point of a single physical element, with: -# (xg, yg) physical Gauss point -# w weight including detJ (or triangle area) -# N length-n_local shape values at the Gauss point -# dN_dx, dN_dy length-n_local physical-coord shape gradients -# -# The same protocol is consumed by assemble_nodal_forces_2d, l2_error_2d, -# and h1_semi_error_2d, so adding a new element type is one rule function. -# --------------------------------------------------------------------------- - -def quad_q1_rule(n_pts=2): - """Element rule for Q1 quads: tensor-product Gauss-Legendre.""" - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"quad_q1_rule: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - def rule(xe, ye): - for xi, wi in zip(xi_pts, w_pts): - for eta, wj in zip(xi_pts, w_pts): - N, dN_dxi, dN_deta = _shape_q1(xi, eta) - J = np.array([[dN_dxi @ xe, dN_dxi @ ye], - [dN_deta @ xe, dN_deta @ ye]]) - detJ = np.linalg.det(J) - Jinv = np.linalg.inv(J) - xg, yg = N @ xe, N @ ye - dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta - dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta - yield xg, yg, wi * wj * detJ, N, dN_dx, dN_dy - return rule - - -# Reference-triangle quadrature, integrating over {(xi,eta): xi,eta>=0, xi+eta<=1}. -# Weights sum to the reference area 1/2; the physical weight is w_ref * 2*area. -_TRI_QUADRATURE = { - 1: (np.array([[1/3, 1/3]]), - np.array([1/2])), - 3: (np.array([[1/6, 1/6], [2/3, 1/6], [1/6, 2/3]]), - np.array([1/6, 1/6, 1/6])), -} - - -def tri_p1_rule(n_pts=1): - """Element rule for P1 triangles: 1-point (centroid) or 3-point Gauss. - - Shape gradients are constant per triangle. Node order in (xe, ye) is the - canonical P1 ordering — corresponds to (N0, N1, N2) = (1-xi-eta, xi, eta). - """ - if n_pts not in _TRI_QUADRATURE: - raise ValueError(f"tri_p1_rule: {n_pts}-point rule not supported") - pts, wts = _TRI_QUADRATURE[n_pts] - - def rule(xe, ye): - x0, x1, x2 = xe - y0, y1, y2 = ye - # Signed double-area: positive for CCW orientation - A2 = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0) - area = abs(A2) / 2.0 - # Constant physical-coord gradients of the P1 shape functions - dN_dx = np.array([(y1 - y2) / A2, (y2 - y0) / A2, (y0 - y1) / A2]) - dN_dy = np.array([(x2 - x1) / A2, (x0 - x2) / A2, (x1 - x0) / A2]) - for (xi, eta), w_ref in zip(pts, wts): - N = np.array([1.0 - xi - eta, xi, eta]) - xg = N[0] * x0 + N[1] * x1 + N[2] * x2 - yg = N[0] * y0 + N[1] * y1 + N[2] * y2 - # w_ref integrates over reference area 1/2; scale by 2*area for physical - yield xg, yg, w_ref * 2.0 * area, N, dN_dx, dN_dy - 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): + """3D Q1 hex shape on [-1, 1]^3. Node order matches SOFA RegularGridTopology: + 0:(-,-,-) 1:(+,-,-) 2:(+,+,-) 3:(-,+,-) + 4:(-,-,+) 5:(+,-,+) 6:(+,+,+) 7:(-,+,+).""" N = 0.125 * np.array([ (1 - xi) * (1 - eta) * (1 - zeta), (1 + xi) * (1 - eta) * (1 - zeta), @@ -167,41 +119,154 @@ def _shape_hex_q1(xi, eta, zeta): return N, dN_dxi, dN_deta, dN_dzeta +# Reference-triangle quadrature, integrating over {(xi,eta): xi,eta>=0, xi+eta<=1}. +# Weights sum to the reference area 1/2; the physical weight is w_ref * 2*area. +_TRI_QUADRATURE = { + 1: (np.array([[1/3, 1/3]]), + np.array([1/2])), + 3: (np.array([[1/6, 1/6], [2/3, 1/6], [1/6, 2/3]]), + np.array([1/6, 1/6, 1/6])), +} + + # --------------------------------------------------------------------------- -# 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) +# Element rules (consume by assemble_nodal_forces, l2_error, h1_semi_error) # --------------------------------------------------------------------------- +def quad_q1_rule(n_pts=2): + """Element rule for Q1 quads: tensor-product Gauss-Legendre on [-1,1]^2.""" + if n_pts not in _GAUSS_LEGENDRE_1D: + raise ValueError(f"quad_q1_rule: {n_pts}-point rule not supported") + xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (4, 2) + x_col, y_col = xe_arr[:, 0], xe_arr[:, 1] + for xi, wi in zip(xi_pts, w_pts): + for eta, wj in zip(xi_pts, w_pts): + N, dN_dxi, dN_deta = _shape_q1(xi, eta) + J = np.array([[dN_dxi @ x_col, dN_dxi @ y_col], + [dN_deta @ x_col, dN_deta @ y_col]]) + detJ = np.linalg.det(J) + Jinv = np.linalg.inv(J) + coords = np.array([N @ x_col, N @ y_col]) + dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta + dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta + dN_phys = np.stack([dN_dx, dN_dy]) # (2, 4) [d, a] + yield coords, wi * wj * detJ, N, dN_phys + return rule + + +def tri_p1_rule(n_pts=1): + """Element rule for P1 triangles: 1-point (centroid) or 3-point Gauss. + + Shape gradients are constant per triangle. Node order in `xe` is the + canonical P1 ordering — corresponds to (N0, N1, N2) = (1-xi-eta, xi, eta). + """ + if n_pts not in _TRI_QUADRATURE: + raise ValueError(f"tri_p1_rule: {n_pts}-point rule not supported") + pts, wts = _TRI_QUADRATURE[n_pts] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (3, 2) + (x0, y0), (x1, y1), (x2, y2) = xe_arr + A2 = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0) + area = abs(A2) / 2.0 + dN_dx = np.array([(y1 - y2) / A2, (y2 - y0) / A2, (y0 - y1) / A2]) + dN_dy = np.array([(x2 - x1) / A2, (x0 - x2) / A2, (x1 - x0) / A2]) + dN_phys = np.stack([dN_dx, dN_dy]) # (2, 3) [d, a] + for (xi, eta), w_ref in zip(pts, wts): + N = np.array([1.0 - xi - eta, xi, eta]) + coords = np.array([N @ xe_arr[:, 0], N @ xe_arr[:, 1]]) + yield coords, w_ref * 2.0 * area, N, dN_phys + return rule + + def hex_q1_rule(n_pts=2): - """Element rule for Q1 hexes: tensor-product Gauss-Legendre over [-1,1]^3.""" + """Element rule for Q1 hexes: tensor-product Gauss-Legendre on [-1,1]^3.""" if n_pts not in _GAUSS_LEGENDRE_1D: raise ValueError(f"hex_q1_rule: {n_pts}-point rule not supported") xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - def rule(xe, ye, ze): + def rule(xe): + xe_arr = np.asarray(xe, float) # (8, 3) + x_col, y_col, z_col = xe_arr[:, 0], xe_arr[:, 1], xe_arr[:, 2] for xi, wi in zip(xi_pts, w_pts): for eta, wj in zip(xi_pts, w_pts): for zeta, wk in zip(xi_pts, w_pts): N, dN_dxi, dN_deta, dN_dzeta = _shape_hex_q1(xi, eta, zeta) J = np.array([ - [dN_dxi @ 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], + [dN_dxi @ x_col, dN_dxi @ y_col, dN_dxi @ z_col], + [dN_deta @ x_col, dN_deta @ y_col, dN_deta @ z_col], + [dN_dzeta @ x_col, dN_dzeta @ y_col, dN_dzeta @ z_col], ]) - detJ = np.linalg.det(J) - Jinv = np.linalg.inv(J) - 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 + detJ = np.linalg.det(J) + Jinv = np.linalg.inv(J) + coords = np.array([N @ x_col, N @ y_col, N @ z_col]) + dN_dx = Jinv[0, 0] * dN_dxi + Jinv[1, 0] * dN_deta + Jinv[2, 0] * dN_dzeta + dN_dy = Jinv[0, 1] * dN_dxi + Jinv[1, 1] * dN_deta + Jinv[2, 1] * dN_dzeta + dN_dz = Jinv[0, 2] * dN_dxi + Jinv[1, 2] * dN_deta + Jinv[2, 2] * dN_dzeta + dN_phys = np.stack([dN_dx, dN_dy, dN_dz]) # (3, 8) [d, a] + yield coords, wi * wj * wk * detJ, N, dN_phys return rule # --------------------------------------------------------------------------- -# FEM assembly +# Facet rules (consumed by assemble_traction) +# --------------------------------------------------------------------------- + +def edge_line_rule(n_pts=2): + """Facet rule for a 2D boundary edge: linear shape, 1D Gauss-Legendre. + + Yields per Gauss point along the edge: (coords (2,), w, N (2,)). + """ + if n_pts not in _GAUSS_LEGENDRE_1D: + raise ValueError(f"edge_line_rule: {n_pts}-point rule not supported") + xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (2, 2) + Le = np.linalg.norm(xe_arr[1] - xe_arr[0]) + for xi, wi in zip(xi_pts, w_pts): + t = 0.5 * (xi + 1.0) + N = np.array([1.0 - t, t]) + coords = N @ xe_arr # (2,) + yield coords, wi * Le / 2.0, N + return rule + + +def quad_face_rule(n_pts=2): + """Facet rule for a 3D boundary quad face: bilinear shape, 2D Gauss. + + Yields per Gauss point on the face: (coords (3,), w, N (4,)). + Quad orientation does not affect the magnitude `dA = ||t_xi × t_eta||`. + """ + if n_pts not in _GAUSS_LEGENDRE_1D: + raise ValueError(f"quad_face_rule: {n_pts}-point rule not supported") + xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] + + def rule(xe): + xe_arr = np.asarray(xe, float) # (4, 3) + for xi, wi in zip(xi_pts, w_pts): + for eta, wj in zip(xi_pts, w_pts): + N = 0.25 * np.array([ + (1 - xi) * (1 - eta), + (1 + xi) * (1 - eta), + (1 + xi) * (1 + eta), + (1 - xi) * (1 + eta), + ]) + dN_dxi = 0.25 * np.array([-(1 - eta), (1 - eta), (1 + eta), -(1 + eta)]) + dN_deta = 0.25 * np.array([-(1 - xi), -(1 + xi), (1 + xi), (1 - xi) ]) + t_xi = dN_dxi @ xe_arr # (3,) + t_eta = dN_deta @ xe_arr + dA = np.linalg.norm(np.cross(t_xi, t_eta)) + coords = N @ xe_arr # (3,) + yield coords, wi * wj * dA, N + return rule + + +# --------------------------------------------------------------------------- +# 1D FEM assembly (scalar-integrand path) # --------------------------------------------------------------------------- def assemble_nodal_forces_1d(f_body, nodes, edges, quadrature): @@ -222,142 +287,68 @@ def assemble_nodal_forces_1d(f_body, nodes, edges, quadrature): return forces -def assemble_traction_2d(traction, nodes, edges, n_pts=2): - """ - Assemble consistent nodal forces from a boundary traction along edges: - - F_a += integral_{edge} N_a(s) * traction(x(s), y(s)) ds - - Linear edge shape functions in physical arc-length. The outward normal is - baked into `traction` by the caller (one lambda per boundary side). - - traction : callable (x, y) -> (Tx, Ty) - nodes : (N, 2) or (N, 3) array — only the first two columns used - edges : iterable of (a, b) node-index pairs along the boundary - n_pts : Gauss points per edge (default 2) - """ - if n_pts not in _GAUSS_LEGENDRE_1D: - raise ValueError(f"assemble_traction_2d: {n_pts}-point rule not supported") - xi_pts, w_pts = _GAUSS_LEGENDRE_1D[n_pts] - - xy = np.asarray(nodes)[:, :2] - F = np.zeros((len(xy), 2)) - for a, b in edges: - x1, y1 = xy[a] - x2, y2 = xy[b] - Le = np.hypot(x2 - x1, y2 - y1) - for xi, wi in zip(xi_pts, w_pts): - t = 0.5 * (xi + 1.0) - xg = (1.0 - t) * x1 + t * x2 - yg = (1.0 - t) * y1 + t * y2 - Tx, Ty = traction(xg, yg) - N0, N1 = 1.0 - t, t - w = wi * Le / 2.0 - F[a, 0] += N0 * Tx * w - F[a, 1] += N0 * Ty * w - F[b, 0] += N1 * Tx * w - F[b, 1] += N1 * Ty * w - return F - +# --------------------------------------------------------------------------- +# Dim-agnostic FEM assembly (2D and 3D) +# --------------------------------------------------------------------------- -def assemble_nodal_forces_2d(f_body, nodes, conn, element_rule): +def assemble_nodal_forces(f_body, nodes, conn, element_rule): """ - Assemble consistent nodal forces on a 2D mesh, agnostic to element type. + Assemble consistent nodal body forces, dim- and element-type agnostic. - F_a = sum_e integral_{Omega_e} N_a(x,y) * f_body(x,y) dx dy + F_a = sum_e integral_{Omega_e} N_a(x) * f_body(x) dx - f_body : callable (x, y) -> (fx, fy) - nodes : (N, 2) or (N, 3) array — only the first two columns used + f_body : callable (*coords) -> dim-tuple of force components + nodes : (N, dim) array conn : iterable of node-index lists, one per element - element_rule : rule from quad_q1_rule(n_pts) / tri_p1_rule(n_pts) + element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule """ - xy = np.asarray(nodes)[:, :2] - F = np.zeros((len(xy), 2)) + nodes = np.asarray(nodes) + dim = nodes.shape[1] + F = np.zeros((len(nodes), dim)) for elem in conn: - xe, ye = xy[elem, 0], xy[elem, 1] - for xg, yg, w, N, _, _ in element_rule(xe, ye): - fx, fy = f_body(xg, yg) - for a, node in enumerate(elem): - F[node, 0] += N[a] * fx * w - F[node, 1] += N[a] * fy * w + # np.asarray forces row-style fancy indexing even when `elem` is a + # Python tuple (which would otherwise trigger multi-axis indexing). + idx = np.asarray(elem) + xe = nodes[idx] # (n_local, dim) + for coords, w, N, _ in element_rule(xe): + f_components = f_body(*coords) + for a, node in enumerate(idx): + for d in range(dim): + F[node, d] += N[a] * f_components[d] * w return F -def assemble_traction_3d(traction, nodes, quads, n_pts=2): +def assemble_traction(traction, nodes, facets, facet_rule): """ - Assemble consistent nodal forces from a boundary traction along quad faces: + Assemble consistent nodal forces from a boundary traction over facets. - F_a += integral_{face} N_a(s,t) * traction(x(s,t), y(s,t), z(s,t)) dA + F_a += integral_{facet} N_a(s) * traction(x(s)) dS - 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`. + The outward unit normal is baked into `traction` by the caller (one + lambda per boundary face); the facet rule only supplies the surface area + element and shape values. - traction : callable (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) + traction : callable (*coords) -> dim-tuple of force components + nodes : (N, dim) array + facets : iterable of node-index tuples (one per boundary facet) + facet_rule : rule from edge_line_rule (2D) / quad_face_rule (3D) """ - 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 + nodes = np.asarray(nodes) + dim = nodes.shape[1] + F = np.zeros((len(nodes), dim)) + for facet in facets: + idx = np.asarray(facet) + xe = nodes[idx] # (n_local, dim) + for coords, w, N in facet_rule(xe): + T = traction(*coords) + for a, node in enumerate(idx): + for d in range(dim): + F[node, d] += N[a] * T[d] * w return F # --------------------------------------------------------------------------- -# Error norms +# 1D error norms # --------------------------------------------------------------------------- def l2_error_1d(nodes, edges, u_h, u_ex, quadrature): @@ -384,126 +375,60 @@ def h1_semi_error_1d(nodes, edges, u_h, du_ex, quadrature): return np.sqrt(total) -def l2_error_2d(nodes, conn, u_h, u_ex, element_rule): - """ - Vector L2 error norm on a 2D mesh, element-type agnostic. - - sqrt( integral || u_h(x,y) - u_ex(x,y) ||^2 dx dy ) over the mesh - - nodes : (N, 2) or (N, 3) array — only the first two columns used - conn : iterable of node-index lists, one per element - u_h : (N, 2) array of nodal displacements (ux, uy) - u_ex : callable (x, y) -> (ux_ex, uy_ex) - element_rule : rule from quad_q1_rule / tri_p1_rule - """ - xy = np.asarray(nodes)[:, :2] - u = np.asarray(u_h) - total = 0.0 - for elem in conn: - xe, ye = xy[elem, 0], xy[elem, 1] - u_loc = u[elem] # (n_local, 2) - for xg, yg, w, N, _, _ in element_rule(xe, ye): - u_h_g = N @ u_loc # (2,) - ux_ex, uy_ex = u_ex(xg, yg) - ex = u_h_g[0] - ux_ex - ey = u_h_g[1] - uy_ex - total += (ex * ex + ey * ey) * w - return np.sqrt(total) - - -def h1_semi_error_2d(nodes, conn, u_h, grad_u_ex, element_rule): - """ - Vector H1 semi-norm error on a 2D mesh, element-type agnostic. - - sqrt( integral || grad u_h - grad u_ex ||_F^2 dx dy ) over the mesh - - nodes : (N, 2) or (N, 3) array — only the first two columns used - conn : iterable of node-index lists, one per element - u_h : (N, 2) array of nodal displacements (ux, uy) - grad_u_ex : callable (x, y) -> 2x2 array - [[dux/dx, dux/dy], [duy/dx, duy/dy]] - element_rule : rule from quad_q1_rule / tri_p1_rule - """ - xy = np.asarray(nodes)[:, :2] - u = np.asarray(u_h) - total = 0.0 - for elem in conn: - xe, ye = xy[elem, 0], xy[elem, 1] - u_loc = u[elem] # (n_local, 2) - for xg, yg, w, _, dN_dx, dN_dy in element_rule(xe, ye): - dux_dx_h = dN_dx @ u_loc[:, 0] - dux_dy_h = dN_dy @ u_loc[:, 0] - duy_dx_h = dN_dx @ u_loc[:, 1] - duy_dy_h = dN_dy @ u_loc[:, 1] - G = grad_u_ex(xg, yg) - total += ( - (dux_dx_h - G[0, 0])**2 + (dux_dy_h - G[0, 1])**2 + - (duy_dx_h - G[1, 0])**2 + (duy_dy_h - G[1, 1])**2 - ) * w - return np.sqrt(total) - +# --------------------------------------------------------------------------- +# Dim-agnostic error norms (2D and 3D) +# --------------------------------------------------------------------------- -def l2_error_3d(nodes, conn, u_h, u_ex, element_rule): +def l2_error(nodes, conn, u_h, u_ex, element_rule): """ - Vector L2 error norm on a 3D mesh, element-type agnostic. + Vector L2 error norm, dim- and element-type agnostic. - sqrt( integral || u_h(x,y,z) - u_ex(x,y,z) ||^2 dx dy dz ) over the mesh + sqrt( integral || u_h(x) - u_ex(x) ||^2 dx ) over the mesh - nodes : (N, 3) array + nodes : (N, dim) 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 + u_h : (N, dim) array of nodal displacements + u_ex : callable (*coords) -> dim-tuple + element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule """ - xyz = np.asarray(nodes)[:, :3] - u = np.asarray(u_h) + nodes = np.asarray(nodes) + 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 + idx = np.asarray(elem) + xe = nodes[idx] # (n_local, dim) + u_loc = u[idx] # (n_local, dim) + for coords, w, N, _ in element_rule(xe): + u_h_g = N @ u_loc # (dim,) + u_ex_g = np.asarray(u_ex(*coords)) # (dim,) + diff = u_h_g - u_ex_g + total += float(np.dot(diff, diff)) * w return np.sqrt(total) -def h1_semi_error_3d(nodes, conn, u_h, grad_u_ex, element_rule): +def h1_semi_error(nodes, conn, u_h, grad_u_ex, element_rule): """ - Vector H1 semi-norm error on a 3D mesh, element-type agnostic. + Vector H1 semi-norm error, dim- and element-type agnostic. - sqrt( integral || grad u_h - grad u_ex ||_F^2 dx dy dz ) over the mesh + sqrt( integral || grad u_h - grad u_ex ||_F^2 dx ) over the mesh - nodes : (N, 3) array + nodes : (N, dim) 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 + u_h : (N, dim) array of nodal displacements + grad_u_ex : callable (*coords) -> (dim, dim) array with [i, j] = du_i/dx_j + element_rule : rule from quad_q1_rule / tri_p1_rule / hex_q1_rule """ - xyz = np.asarray(nodes)[:, :3] - u = np.asarray(u_h) + nodes = np.asarray(nodes) + 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 + idx = np.asarray(elem) + xe = nodes[idx] # (n_local, dim) + u_loc = u[idx] # (n_local, dim) + for coords, w, _, dN_phys in element_rule(xe): + # grad_uh[i, d] = sum_a u_i[a] * dN_a/dx_d + grad_uh = (dN_phys @ u_loc).T # (dim, dim) [i, d] + grad_ue = np.asarray(grad_u_ex(*coords)) # [i, j] = du_i/dx_j + diff = grad_uh - grad_ue + total += float(np.sum(diff * diff)) * w return np.sqrt(total) diff --git a/examples/Freefem/MMS/mms_case.py b/examples/Freefem/MMS/mms_case.py new file mode 100644 index 00000000..8f426ace --- /dev/null +++ b/examples/Freefem/MMS/mms_case.py @@ -0,0 +1,15 @@ +"""Cross-dim base class for manufactured-solution cases. + +Dim-specific bases (`MMSCase1D`, `MMSCase2D`, `MMSCase3D`) inherit from +`MMSCase` and declare the abstract solution/derivative/source/BC methods +plus their `source_quadrature_*` placeholders. The case files +(`cubic.py`, `trigonometric.py`, `sinus_neumann.py`, …) inherit from the +dim-specific bases and override the abstract methods explicitly. +""" + +from abc import ABC + + +class MMSCase(ABC): + name = None # case identifier (must match the params.json key) + plot_label = None # LaTeX label for the exact solution diff --git a/examples/Freefem/MMS/output.py b/examples/Freefem/MMS/output.py new file mode 100644 index 00000000..29d7c34c --- /dev/null +++ b/examples/Freefem/MMS/output.py @@ -0,0 +1,105 @@ +"""Shared output helpers across MMS 1D/2D/3D drivers (tables + plots).""" + +import os +import numpy as np +import matplotlib.pyplot as plt + +from plot_utils import annotate_convergence_rates + + +_AXES = ("x", "y", "z") + + +def write_solution_table(stem, coords, u_h, u_ex, results_dir, error_dict): + """ + Write per-node solution table and error summary to /.txt. + + coords : (N,) or (N, dim) array of node coordinates + u_h : (N,) or (N, dim) array of nodal displacements + u_ex : callable (*coords_row) -> scalar (dim=1) or dim-tuple + error_dict : ordered dict of {label: value} for error summary lines + + Column count adjusts with dim. For 1D the displacement column header is + `ux_h` and the error header is `err_x` (uniform with 2D/3D). + """ + coords = np.asarray(coords) + if coords.ndim == 1: + coords = coords.reshape(-1, 1) + u_h = np.asarray(u_h) + if u_h.ndim == 1: + u_h = u_h.reshape(-1, 1) + dim = coords.shape[1] + assert u_h.shape[1] == dim, "u_h column count must match coord dim" + axes = _AXES[:dim] + + coord_hdr = " | ".join(f"{ax:>10}" for ax in axes) + disp_hdr = " | ".join(f"{'u'+ax+'_h':>15}" for ax in axes) + exact_hdr = " | ".join(f"{'u'+ax+'_ex':>15}" for ax in axes) + err_hdr = " | ".join(f"{'err_'+ax:>15}" for ax in axes) + header = f"{coord_hdr} | {disp_hdr} | {exact_hdr} | {err_hdr}" + + os.makedirs(results_dir, exist_ok=True) + path = os.path.join(results_dir, f"{stem}.txt") + with open(path, "w") as f: + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for c_row, uh_row in zip(coords, u_h): + ue = u_ex(*c_row) + ue = (ue,) if np.isscalar(ue) else tuple(ue) + coord_str = " | ".join(f"{x:10.4f}" for x in c_row) + disp_str = " | ".join(f"{x:15.6e}" for x in uh_row) + exact_str = " | ".join(f"{x:15.6e}" for x in ue) + err_str = " | ".join(f"{abs(uh_row[k] - ue[k]):15.6e}" + for k in range(dim)) + f.write(f"{coord_str} | {disp_str} | {exact_str} | {err_str}\n") + f.write("\n") + for label, val in error_dict.items(): + f.write(f"{label:12s} = {val:.6e}\n") + + +def write_convergence_table(stem, rows, results_dir): + """ + Write convergence table to /.txt. + + rows : list of dicts with keys 'nx', 'h', and one key per error column. + Rate columns are strings (empty for the first row). + """ + os.makedirs(results_dir, exist_ok=True) + path = os.path.join(results_dir, f"{stem}.txt") + err_keys = [k for k in rows[0] if k not in ("nx", "h")] + header = f"{'nx':>6} | {'h':>10}" + "".join(f" | {k:>16}" for k in err_keys) + with open(path, "w") as f: + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for row in rows: + line = f"{row['nx']:6d} | {row['h']:10.4f}" + for k in err_keys: + v = row[k] + line += f" | {v:16.6e}" if isinstance(v, float) else f" | {v:>16}" + f.write(line + "\n") + + +def plot_convergence(stem, hs, series, title, results_dir, ylabel="Error"): + """ + Save log-log convergence plot to /.png. + + series : list of {"label", "errors", "style"?} dicts + Per-segment convergence rates are annotated above each line segment. + """ + os.makedirs(results_dir, exist_ok=True) + h_arr = np.array(hs) + default = ["bo-", "rs--", "g^:", "m^-"] + fig, ax = plt.subplots(figsize=(8, 5)) + for i, s in enumerate(series): + style = s.get("style", default[i % len(default)]) + e_arr = np.array(s["errors"]) + ax.loglog(h_arr, e_arr, style, label=s["label"], linewidth=2, markersize=7) + annotate_convergence_rates(ax, h_arr, e_arr) + ax.set_xlabel("h") + ax.set_ylabel(ylabel) + ax.set_title(title) + ax.legend() + ax.grid(True, alpha=0.3, which="both") + fig.tight_layout() + fig.savefig(os.path.join(results_dir, f"{stem}.png"), dpi=150) + plt.close(fig) diff --git a/examples/Freefem/MMS/scene.py b/examples/Freefem/MMS/scene.py new file mode 100644 index 00000000..bcd390fd --- /dev/null +++ b/examples/Freefem/MMS/scene.py @@ -0,0 +1,40 @@ +"""Shared SOFA scene helpers (controllers) across MMS 1D/2D/3D drivers.""" + +import Sofa +import Sofa.Core + + +class NodalForceAssembler(Sofa.Core.Controller): + """Fill a placeholder ConstantForceField after Sofa init has run. + + SOFA topology components (RegularGridTopology, + EdgeSetTopologyContainer, QuadSetTopologyContainer, + HexahedronSetTopologyContainer, …) are only populated after init runs, + not during Python scene-build time. This controller defers nodal-force + assembly to onSimulationInitDoneEvent, where it reads rest positions + off the MechanicalObject and hands them to `compute_forces` together + with the topology container. The returned array is written into the + force field. + + Parameters + ---------- + dofs : SOFA MechanicalObject with rest_position + topology : SOFA topology container (read by compute_forces) + force_field : SOFA ConstantForceField to populate + compute_forces : callable (nodes, topology) -> ndarray sized like + force_field.forces + """ + + def __init__(self, dofs, topology, force_field, compute_forces, + *args, **kwargs): + super().__init__(*args, **kwargs) + self.dofs = dofs + self.topology = topology + self.force_field = force_field + self.compute_forces = compute_forces + + def onSimulationInitDoneEvent(self, event): + nodes = self.dofs.rest_position.array().copy() + F = self.compute_forces(nodes, self.topology) + with self.force_field.forces.writeableArray() as forces: + forces[:] = F