From 5f08ced3fdf4582fc5b17fb86f13e9c26313f523 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 25 May 2026 20:41:01 +0200 Subject: [PATCH 1/4] Refactor 2D and revise 1D - Unify fem.py API for 2D - Factorize the Neumann BC traction computation - ABC introduction - Refactored cubic and trigonometric - Use SOFA topology - Restructured params.json --- examples/Freefem/MMS/1D/bar.py | 26 +- examples/Freefem/MMS/1D/params.json | 16 +- examples/Freefem/MMS/1D/run_convergence.py | 80 +- examples/Freefem/MMS/2D/beam.py | 490 ++++++++++ examples/Freefem/MMS/2D/beam_solution.py | 13 + examples/Freefem/MMS/2D/cubic.py | 56 ++ .../Freefem/MMS/2D/manufactured_solution.py | 56 ++ examples/Freefem/MMS/2D/mms_2d.py | 667 -------------- examples/Freefem/MMS/2D/params.json | 17 + examples/Freefem/MMS/2D/run_convergence.py | 163 ++++ .../Freefem/MMS/2D/trigonometric-cubic-mms.py | 836 ------------------ examples/Freefem/MMS/2D/trigonometric.py | 62 ++ examples/Freefem/MMS/fem.py | 224 ++++- examples/Freefem/MMS/plot_utils.py | 16 + 14 files changed, 1163 insertions(+), 1559 deletions(-) create mode 100644 examples/Freefem/MMS/2D/beam.py create mode 100644 examples/Freefem/MMS/2D/beam_solution.py create mode 100644 examples/Freefem/MMS/2D/cubic.py create mode 100644 examples/Freefem/MMS/2D/manufactured_solution.py delete mode 100644 examples/Freefem/MMS/2D/mms_2d.py create mode 100644 examples/Freefem/MMS/2D/params.json create mode 100644 examples/Freefem/MMS/2D/run_convergence.py delete mode 100644 examples/Freefem/MMS/2D/trigonometric-cubic-mms.py create mode 100644 examples/Freefem/MMS/2D/trigonometric.py create mode 100644 examples/Freefem/MMS/plot_utils.py diff --git a/examples/Freefem/MMS/1D/bar.py b/examples/Freefem/MMS/1D/bar.py index 57493755..8a6c2f41 100644 --- a/examples/Freefem/MMS/1D/bar.py +++ b/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", @@ -141,7 +141,7 @@ 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 +211,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) diff --git a/examples/Freefem/MMS/1D/params.json b/examples/Freefem/MMS/1D/params.json index 4c846109..c094267f 100644 --- a/examples/Freefem/MMS/1D/params.json +++ b/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] + } } } diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py index 70f64bc6..b68615af 100644 --- a/examples/Freefem/MMS/1D/run_convergence.py +++ b/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/_convergence.txt. + 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"{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,49 @@ 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", ref_lines=None): """ - Save log-log convergence plot to results/_convergence.png. + Save log-log convergence plot to results/.png. - error_series : dict { label: array-like of errors } + series : list of {"label", "errors", "style"?} dicts + ref_lines : optional list of {"slope", "anchor", "label"?} dicts; + "anchor" names the series whose first point anchors the line. 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) + anchors = {} + 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) + anchors[s["label"]] = e_arr + for rl in (ref_lines or []): + anchor = anchors.get(rl["anchor"]) + if anchor is None: + continue + slope = rl["slope"] + label = rl.get("label", f"O(h^{slope})") + h_ref = np.array([h_arr[0], h_arr[-1]]) + ax.loglog(h_ref, anchor[0] * (h_ref / h_arr[0])**slope, ":", + color="gray", lw=1.2, label=label) 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)}) diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py new file mode 100644 index 00000000..0e2046db --- /dev/null +++ b/examples/Freefem/MMS/2D/beam.py @@ -0,0 +1,490 @@ +import json +import numpy as np +import matplotlib.pyplot as plt +import os +import sys + +import Sofa +import Sofa.Core +import Sofa.Simulation +import SofaRuntime + +# Make the parent MMS/ directory importable so we can pull in fem.py. +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +from fem import ( + assemble_nodal_forces_2d, + assemble_traction_2d, + l2_error_2d, + h1_semi_error_2d, + quad_q1_rule, + tri_p1_rule, +) +from beam_solution import BeamSolution2D + +RESULTS_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results") + + +# --------------------------------------------------------------------------- +# Parameters +# --------------------------------------------------------------------------- + +def load_params(path=None): + if path is None: + path = os.path.join(os.path.dirname(os.path.abspath(__file__)), + "params.json") + with open(path) as f: + return json.load(f) + +# ============== Mesh helpers ============================ + + +def get_nodes_2d(L, nx, ny, dim="2d"): + dx = L / (nx - 1) + dy = L / (ny - 1) + if dim == "3d": + pts = [[i*dx, j*dy, 0.0] for j in range(ny) for i in range(nx)] + else: + pts = [[i*dx, j*dy] for j in range(ny) for i in range(nx)] + return np.array(pts, dtype=float) + + +def _dim_template(dim): + return "Vec3d" if dim == "3d" else "Vec2d" + + +def _add_dirichlet(Beam, nodes_2d, L, dim, with_visual): + """Partial Dirichlet: ux=0 on x-faces, uy=0 on y-faces (and z=0 in 3d).""" + tmpl = _dim_template(dim) + xy = nodes_2d[:, :2] + eps = 1e-10 + + idx_ux0 = [k for k, (xk, _) in enumerate(xy) if xk < eps or xk > L - eps] + idx_uy0 = [k for k, (_, yk) in enumerate(xy) if yk < eps or yk > L - eps] + + if dim == "2d": + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1") + else: + all_idx = " ".join(map(str, range(len(nodes_2d)))) + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_z", template=tmpl, + indices=all_idx, + fixedDirections="0 0 1") + + +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.ELEMENT_RULE) + + 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)) + + @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)) + + @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`. +# --------------------------------------------------------------------------- + +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 + + +# --------------------------------------------------------------------------- +# SOFA scene +# --------------------------------------------------------------------------- + +def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, + nx=10, ny=10, with_visual=True, dim="2d"): + """Build a SOFA scene for `mms` on the 2D `element` strategy. + + dim : "2d" → Vec2d template / plane stress + "3d" → Vec3d template / plane strain (z coordinate fixed at 0) + Returns (dofs, topology). Nodes and connectivity become available + after `Sofa.Simulation.init(root)` runs, via + `dofs.rest_position.array()` and `element.read_connectivity(topology)`. + """ + tmpl = _dim_template(dim) + + rootNode.addObject("RequiredPlugin", pluginName=[ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Grid", + "Sofa.Component.Topology.Container.Dynamic", + "Sofa.Component.Topology.Mapping", + "Sofa.Component.Visual", + ]) + rootNode.addObject("DefaultAnimationLoop") + if with_visual: + rootNode.addObject("VisualStyle", + displayFlags="showBehaviorModels showForceFields") + + # Nodes still computed in Python so MechanicalObject's Vec2d/Vec3d template + # receives the right per-point dimension; SOFA-side positions on the + # topology container (via @../Grid/grid.position) are kept in sync by the + # regular-grid construction. + nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) + + # Grid must be added *before* Beam: SOFA inits children in insertion + # order, and Beam's topology container resolves `@../Grid/grid.position` + # during its own init — so Grid has to be initialised first. + Grid = rootNode.addChild("Grid") + Grid.addObject("RegularGridTopology", name="grid", + nx=nx, ny=ny, nz=1, + min=[0.0, 0.0, 0.0], max=[L, L, 0.0]) + + Beam = rootNode.addChild("Beam") + Beam.addObject("StaticSolver", name="staticSolver", printLog=False) + Beam.addObject("NewtonRaphsonSolver", name="newtonSolver", + maxNbIterationsNewton=1, + absoluteResidualStoppingThreshold=1e-10, + printLog=False) + Beam.addObject("SparseLDLSolver", name="linearSolver", + template="CompressedRowSparseMatrixd") + + dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, + position=nodes_2d.tolist(), + showObject=with_visual, showObjectScale=0.005 * L) + + topology = element.add_topology(Beam) + + Beam.addObject("LinearSmallStrainFEMForceField", name="FEM", template=tmpl, + youngModulus=E, poissonRatio=nu, topology="@topology") + + _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) + + # Placeholder force field filled in by the controller after init. + n_nodes = len(nodes_2d) + n_comp = 3 if dim == "3d" else 2 + force_field = Beam.addObject("ConstantForceField", name="MMS_forces", + template=tmpl, + indices=list(range(n_nodes)), + forces=[[0.0] * n_comp] * n_nodes) + + 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, + name="nodalForceAssembler")) + + return dofs, topology + + +# ───────────────────────────────────────────────────────────────────────────── +# Simulation runner +# ───────────────────────────────────────────────────────────────────────────── + +def solve_beam(elem, mms, L, E, nu, nx, ny, dim="2d"): + """Build, init, and run one static step. Returns a BeamSolution2D snapshot.""" + root = Sofa.Core.Node("root") + dofs, topology = build_beam_scene( + root, mms, elem, L=L, E=E, nu=nu, nx=nx, ny=ny, with_visual=False, dim=dim + ) + Sofa.Simulation.init(root) + # Read topology back from SOFA now that init has populated it. + nodes_2d = dofs.rest_position.array().copy() + conn = elem.read_connectivity(topology) + pos0 = dofs.position.array().copy() + Sofa.Simulation.animate(root, root.dt.value) + pos1 = dofs.position.array().copy() + Sofa.Simulation.unload(root) + # Extract x and y displacement columns (works for both Vec2d and Vec3d) + ux = pos1[:, 0] - pos0[:, 0] + uy = pos1[:, 1] - pos0[:, 1] + return BeamSolution2D(nodes=nodes_2d, conn=conn, ux=ux, uy=uy) + + +# --------------------------------------------------------------------------- +# 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) + xy = sol.nodes[:, :2] + mid_j = (ny - 1) // 2 + sl = slice(mid_j * nx, mid_j * nx + nx) + yc = xy[mid_j * nx, 1] + xf = np.linspace(0, L, 300) + + fig, axes = plt.subplots(1, 2, figsize=(12, 4)) + for ax, u_sofa, u_fn, lbl, fmt in zip( + axes, + [sol.ux[sl], sol.uy[sl]], + [lambda x: mms.u_ex(x, yc, L)[0], lambda x: mms.u_ex(x, yc, L)[1]], + [r"$u_x$", r"$u_y$"], + ["o-", "s-"], + ): + ax.plot(xy[sl, 0], u_sofa, fmt, color="tab:green", + label=f"SOFA {label} [{dim}]", ms=5) + ax.plot(xf, u_fn(xf), "--", color="tab:blue", label="MMS exact") + ax.set_xlabel("x"); ax.set_ylabel(lbl) + ax.legend(); ax.grid(True, alpha=0.3) + fig.suptitle(f"{mms.name} — {label} [{dim} / {hyp}] " + f"nu={nu} nx={nx} |L2={l2:.2e} H1={h1:.2e}") + fig.tight_layout() + fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150) + plt.close(fig) + + +def plot_solution_fields(stem, sol, elem, mms, L, nx, label, dim, hyp, nu): + """Save 2-D field colour maps (ux, uy, error) to results/.png.""" + os.makedirs(RESULTS_DIR, exist_ok=True) + xy = sol.nodes[:, :2] + ux_ref, uy_ref = mms.u_ex(xy[:, 0], xy[:, 1], L) + tris_plot = elem.to_triangles(sol.conn) + x, y = xy[:, 0], xy[:, 1] + + fig, axes = plt.subplots(2, 3, figsize=(15, 8)) + for ax, data, title, cmap in [ + (axes[0,0], sol.ux, r"$u_x$ SOFA", "gray"), + (axes[0,1], ux_ref, r"$u_x$ MMS", "gray"), + (axes[0,2], np.abs(sol.ux-ux_ref), r"$|u_x - u_x^{MMS}|$", "gray_r"), + (axes[1,0], sol.uy, r"$u_y$ SOFA", "gray"), + (axes[1,1], uy_ref, r"$u_y$ MMS", "gray"), + (axes[1,2], np.abs(sol.uy-uy_ref), r"$|u_y - u_y^{MMS}|$", "gray_r"), + ]: + tc = ax.tricontourf(x, y, tris_plot.tolist(), data, levels=20, cmap=cmap) + ax.triplot(x, y, tris_plot.tolist(), "k-", lw=0.3, alpha=0.4) + plt.colorbar(tc, ax=ax, shrink=0.8) + ax.set_title(title); ax.set_aspect("equal") + ax.set_xlabel("x"); ax.set_ylabel("y") + fig.suptitle(f"Fields 2D — {label} [{dim} / {hyp}] " + f"{mms.name} nu={nu} nx={nx}") + fig.tight_layout() + fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150) + plt.close(fig) + + +# ───────────────────────────────────────────────────────────────────────────── +# Single-case driver +# ───────────────────────────────────────────────────────────────────────────── + +def run_reference_scene(elem, mms): + """Solve one MMS case at the reference mesh, write the solution table and plot. + + All parameters come from params.json (top-level + `reference` block). + """ + cfg = load_params() + ref = cfg["reference"] + L, E = cfg["length"], cfg["youngModulus"] + nu, dim = ref["nu"], ref["dim"] + nx = ny = ref["nx"] + hyp = "plane strain" if dim == "3d" else "plane stress" + + 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) + + label = elem.LABEL + tag = label.replace(" ", "_") + 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, + lambda xi, yi: mms.u_ex(xi, yi, L), + {"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, + label, dim, hyp, nu) + + +def case_scene(mms, element): + """Return a `createScene(rootNode)` bound to this MMS and element type. + + All parameters come from params.json (top-level + `reference` block). Each + case file exposes: + createScene = case_scene(mms, element_quad) + so that `runSofa cubic.py` loads the default scene. + """ + def createScene(rootNode): + cfg = load_params() + ref = cfg["reference"] + build_beam_scene(rootNode, mms, element, + L=cfg["length"], E=cfg["youngModulus"], + nu=ref["nu"], nx=ref["nx"], ny=ref["nx"], + 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/beam_solution.py b/examples/Freefem/MMS/2D/beam_solution.py new file mode 100644 index 00000000..4681afed --- /dev/null +++ b/examples/Freefem/MMS/2D/beam_solution.py @@ -0,0 +1,13 @@ +"""Snapshot of one 2D solve: mesh + connectivity + displacement field.""" + +from dataclasses import dataclass + +import numpy as np + + +@dataclass +class BeamSolution2D: + nodes : np.ndarray # (N, 2) or (N, 3) + conn : np.ndarray # element connectivity (element-type-specific) + ux : np.ndarray # (N,) x-displacement + uy : np.ndarray # (N,) y-displacement diff --git a/examples/Freefem/MMS/2D/cubic.py b/examples/Freefem/MMS/2D/cubic.py new file mode 100644 index 00000000..43c17afc --- /dev/null +++ b/examples/Freefem/MMS/2D/cubic.py @@ -0,0 +1,56 @@ +""" +Cubic 2D MMS on [0,L]^2 with linear-elasticity constitutive law: + + u_ex(x, y) = ( x (L - x)(x + y), + y (L - y)(y - x) ) + + sigma = lambda tr(eps) I + 2 mu eps, + with (lambda, mu) selected per dim (plane stress vs plane strain). +""" + +import numpy as np + +from manufactured_solution import MMSCase2D, lame + + +class Cubic(MMSCase2D): + name = "cubic" + plot_label = r"$u_x = x(L-x)(x+y),\ u_y = y(L-y)(y-x)$" + + def u_ex(self, x, y, L): + return (x * (L - x) * (x + y), + y * (L - y) * (y - x)) + + def grad_u_ex(self, x, y, L): + dux_dx = (L - 2*x) * (x + y) + x * (L - x) + dux_dy = x * (L - x) + duy_dx = -y * (L - y) + duy_dy = (L - 2*y) * (y - x) + y * (L - y) + return np.array([[dux_dx, dux_dy], + [duy_dx, duy_dy]]) + + def source(self, x, y, E, nu, L, dim): + lam, mu = lame(E, nu, dim) + d2ux_dxx = 2*L - 6*x - 2*y + d2ux_dxy = L - 2*x + d2ux_dyy = np.zeros_like(np.asarray(x, float)) + d2uy_dxx = np.zeros_like(np.asarray(x, float)) + d2uy_dxy = -(L - 2*y) + d2uy_dyy = 2*L - 6*y + 2*x + fx = -((lam + 2*mu) * d2ux_dxx + lam * d2uy_dxy + + mu * (d2ux_dyy + d2uy_dxy)) + fy = -(mu * (d2ux_dxy + d2uy_dxx) + lam * d2ux_dxy + + (lam + 2*mu) * d2uy_dyy) + return (fx, fy) + + +from beam import case_scene, element_quad, element_tri + +mms = Cubic() +createScene = case_scene(mms, element_quad) + + +if __name__ == "__main__": + from beam import run_reference_scene + + run_reference_scene(element_quad, mms) diff --git a/examples/Freefem/MMS/2D/manufactured_solution.py b/examples/Freefem/MMS/2D/manufactured_solution.py new file mode 100644 index 00000000..88ce9408 --- /dev/null +++ b/examples/Freefem/MMS/2D/manufactured_solution.py @@ -0,0 +1,56 @@ +"""Abstract base class for 2D MMS cases (Cartesian domain [0,L]^2).""" + +from abc import ABC, abstractmethod + + +def lame(E, nu, dim): + """Lamé parameters (lambda, mu) for linear elasticity. + + dim="2d" → plane stress: lambda = E nu / (1 - nu^2) + dim="3d" → plane strain: lambda = E nu / ((1 + nu)(1 - 2 nu)) + mu = E / (2 (1 + nu)) (same in both branches) + """ + if dim == "2d": + lam = E * nu / (1.0 - nu**2) + else: + lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) + mu = E / (2.0 * (1.0 + nu)) + 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 + + @abstractmethod + def u_ex(self, x, y, L): + """Exact solution: returns (ux, uy).""" + + @abstractmethod + def grad_u_ex(self, x, y, L): + """Exact displacement gradient: returns 2x2 array + [[dux/dx, dux/dy], [duy/dx, duy/dy]].""" + + @abstractmethod + def source(self, x, y, E, nu, L, dim): + """Body force: returns (fx, fy). + + `dim` selects the constitutive branch: + "2d" → plane stress, "3d" → plane strain. + """ + + def traction(self, x, y, nx, ny, E, nu, L, dim): + """Traction on a face with outward unit normal (nx, ny): returns (Tx, Ty). + + Derived from `grad_u_ex` via sigma = lambda tr(eps) I + 2 mu eps and + T = sigma . n. `dim` selects plane stress / plane strain (see `source`). + """ + lam, mu = lame(E, nu, dim) + G = self.grad_u_ex(x, y, L) + dux_dx, dux_dy = G[0, 0], G[0, 1] + duy_dx, duy_dy = G[1, 0], G[1, 1] + sxx = (lam + 2.0 * mu) * dux_dx + lam * duy_dy + syy = lam * dux_dx + (lam + 2.0 * mu) * duy_dy + sxy = mu * (dux_dy + duy_dx) + return (sxx * nx + sxy * ny, + sxy * nx + syy * ny) diff --git a/examples/Freefem/MMS/2D/mms_2d.py b/examples/Freefem/MMS/2D/mms_2d.py deleted file mode 100644 index 75d5796d..00000000 --- a/examples/Freefem/MMS/2D/mms_2d.py +++ /dev/null @@ -1,667 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import os - -import Sofa -import Sofa.Core -import Sofa.Simulation -import SofaRuntime - -RESULTS_DIR = "results_mms2d_Q1P1" -os.makedirs(RESULTS_DIR, exist_ok=True) - -GAUSS_PTS = np.array([-1.0 / np.sqrt(3), 1.0 / np.sqrt(3)]) -GAUSS_WTS = np.array([1.0, 1.0]) - - - -#============ MMS manufactured solution & derivatives ===================== - - -def ux_mms(x, y, L): return x**2 * (L - x) / L**2 - -def uy_mms(x, y, L): return x * (L - x) * y / L**2 - -def dux_dx(x, y, L): return (2*x*L - 3*x**2) / L**2 - -def dux_dy(x, y, L): return np.zeros_like(np.asarray(x, float)) - -def duy_dx(x, y, L): return (L - 2*x) * y / L**2 - -def duy_dy(x, y, L): return x * (L - x) / L**2 - -def sigma_xx(x, y, E, L): return E * dux_dx(x, y, L) - -def sigma_yy(x, y, E, L): return E * duy_dy(x, y, L) - -def sigma_xy(x, y, E, L): return E * 0.5 * duy_dx(x, y, L) - -def fx_body(x, y, E, L): return -(E*(2*L - 6*x)/L**2 + E*(L - 2*x)/(2*L**2)) - -def fy_body(x, y, E, L): return E * y / L**2 - -def traction(x, y, nx_c, ny_c, E, L): - sxx = sigma_xx(x, y, E, L) - syy = sigma_yy(x, y, E, L) - sxy = sigma_xy(x, y, E, L) - return sxx*nx_c + sxy*ny_c, sxy*nx_c + syy*ny_c - - - -# ============== Mesh helpers ============================ - - -def get_nodes_2d(L, nx, ny, dim="2d"): - dx = L / (nx - 1) - dy = L / (ny - 1) - if dim == "3d": - pts = [[i*dx, j*dy, 0.0] for j in range(ny) for i in range(nx)] - else: - pts = [[i*dx, j*dy] for j in range(ny) for i in range(nx)] - return np.array(pts, dtype=float) - - -def _dim_template(dim): - - return "Vec3d" if dim == "3d" else "Vec2d" - - - -# =========== Q1 quad elements ========================= - - -class _QuadElement: - LABEL = "Q1 quad" - - @staticmethod - def get_quads(nx, ny): - quads = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - quads.append([k00, k10, k11, k01]) - return np.array(quads) - - @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) - - @staticmethod - def _shape_q1(xi, eta): - 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)]) - return N, dN_dxi, dN_deta - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nx, ny): - """Always uses the first two coordinate columns (x, y).""" - quads = _QuadElement.get_quads(nx, ny) - # Work with xy columns regardless of 2d/3d storage - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - - # Body forces -- Gauss quadrature over each quad - for quad in quads: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape_q1(xi, eta) - J = np.array([[dN_dxi@xe, dN_dxi@ye], - [dN_deta@xe, dN_deta@ye]]) - detJ = np.linalg.det(J) - xg, yg = N@xe, N@ye - fx, fy = fx_body(xg, yg, E, L), fy_body(xg, yg, E, L) - w = wi * wj * detJ - for a, node in enumerate(quad): - F[node, 0] += N[a] * fx * w - F[node, 1] += N[a] * fy * w - - # Neumann BC — bottom edge (y=0) - for i in range(nx - 1): - k0, k1 = i, i + 1 - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5*(xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, 0., -1., E, L) - N0, N1 = 0.5*(1-xi), 0.5*(1+xi) - w = wi * Le / 2.0 - F[k0, 0] += N0*Tx*w; F[k0, 1] += N0*Ty*w - F[k1, 0] += N1*Tx*w; F[k1, 1] += N1*Ty*w - - # Neumann BC — top edge (y=L) - for i in range(nx - 1): - k0 = (ny-1)*nx + i - k1 = (ny-1)*nx + i + 1 - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5*(xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, 0., +1., E, L) - N0, N1 = 0.5*(1-xi), 0.5*(1+xi) - w = wi * Le / 2.0 - F[k0, 0] += N0*Tx*w; F[k0, 1] += N0*Ty*w - F[k1, 0] += N1*Tx*w; F[k1, 1] += N1*Ty*w - return F # shape (N, 2) — z-component added in create_sofa_scene if needed - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nx=10, ny=10, - with_visual=True, dim="2d"): - """ - Build SOFA scene for Q1 quads. - dim : "2d" → Vec2d template - "3d" → Vec3d template, z coordinate fixed at 0 - Returns (dofs, nodes_2d, quads). - """ - tmpl = _dim_template(dim) - - rootNode.addObject("RequiredPlugin", pluginName="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - quads_np = _QuadElement.get_quads(nx, ny) - - Beam = rootNode.addChild("Beam") - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject('NewtonRaphsonSolver', - name="newtonSolver", - maxNbIterationsNewton=1, - absoluteResidualStoppingThreshold=1e-10, - printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), showObject=with_visual, - showObjectScale=0.005*L) - - Beam.addObject("QuadSetTopologyContainer", name="topology", - quads=quads_np.tolist()) - Beam.addObject("QuadSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", name="FEM", template=tmpl, - youngModulus=E, poissonRatio=0.0, topology="@topology") - - e = 1e-4 * L / max(nx-1, ny-1) - if dim == "3d": - box_left = [-e, -e, -e, e, L+e, e] - box_right = [L-e, -e, -e, L+e, L+e, e] - else: - box_left = [-e, -e, 0, e, L+e, 0] - box_right = [L-e, -e, 0, L+e, L+e, 0] - - Beam.addObject("BoxROI", name="box_left", template=tmpl, - box=box_left, drawBoxes=with_visual) - Beam.addObject("FixedProjectiveConstraint", name="fix_left", template=tmpl, - indices="@box_left.indices") - Beam.addObject("BoxROI", name="box_right", template=tmpl, - box=box_right, drawBoxes=with_visual) - Beam.addObject("FixedProjectiveConstraint", name="fix_right", template=tmpl, - indices="@box_right.indices") - - # Build nodal force vector; pad with z=0 column for Vec3d - F_xy = _QuadElement.compute_nodal_forces(nodes_2d, L, E, nx, ny) - if dim == "3d": - F_all = np.hstack([F_xy, np.zeros((len(F_xy), 1))]) - else: - F_all = F_xy - - idx = " ".join(str(k) for k in range(len(nodes_2d))) - if dim == "3d": - frc = " ".join(f"{F_all[k,0]} {F_all[k,1]} {F_all[k,2]}" - for k in range(len(nodes_2d))) - else: - frc = " ".join(f"{F_all[k,0]} {F_all[k,1]}" - for k in range(len(nodes_2d))) - - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - return dofs, nodes_2d, quads_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, quads_np): - xy = nodes_2d[:, :2] - err2 = 0.0 - for quad in quads_np: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape_q1(xi, eta) - J = np.array([[dN_dxi@xe, dN_dxi@ye], - [dN_deta@xe, dN_deta@ye]]) - detJ = np.linalg.det(J) - xg, yg = N@xe, N@ye - ux_h = N @ ux[quad]; uy_h = N @ uy[quad] - ex = ux_h - ux_mms(xg, yg, L) - ey = uy_h - uy_mms(xg, yg, L) - err2 += (ex**2 + ey**2) * wi * wj * detJ - return np.sqrt(err2) - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, quads_np, **_ignored): - xy = nodes_2d[:, :2] - err2 = 0.0 - for quad in quads_np: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._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 - dux_dx_h = dN_dx @ ux[quad]; dux_dy_h = dN_dy @ ux[quad] - duy_dx_h = dN_dx @ uy[quad]; duy_dy_h = dN_dy @ uy[quad] - err2 += ( - (dux_dx_h - dux_dx(xg,yg,L))**2 + - (dux_dy_h - dux_dy(xg,yg,L))**2 + - (duy_dx_h - duy_dx(xg,yg,L))**2 + - (duy_dy_h - duy_dy(xg,yg,L))**2 - ) * wi * wj * detJ - return np.sqrt(err2) - - -# P1 triangle elements - - -class _TriElement: - LABEL = "P1 tri" - - @staticmethod - def get_triangles(nx, ny): - tris = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - tris += [[k00, k10, k11], [k00, k11, k01]] - return np.array(tris) - - @staticmethod - def to_triangles(conn): - return conn - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nx, ny): - """Always uses the first two coordinate columns (x, y).""" - dx, dy = L / (nx - 1), L / (ny - 1) - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - eps = 1e-10 - - # Body forces — trapezoidal nodal quadrature - for k, (xk, yk) in enumerate(xy): - i = round(xk / dx); j = round(yk / dy) - wx = dx/2 if (i == 0 or i == nx-1) else dx - wy = dy/2 if (j == 0 or j == ny-1) else dy - F[k, 0] += fx_body(xk, yk, E, L) * wx * wy - F[k, 1] += fy_body(xk, yk, E, L) * wx * wy - - # Neumann BC — bottom edge (y=0) - for i in range(nx): - k = i - xk, yk = xy[k] - if xk < eps or xk > L - eps: continue - wx = dx/2 if (i == 0 or i == nx-1) else dx - Tx, Ty = traction(xk, yk, 0., -1., E, L) - F[k, 0] += Tx * wx; F[k, 1] += Ty * wx - - # Neumann BC — top edge (y=L) - for i in range(nx): - k = (ny-1)*nx + i - xk, yk = xy[k] - if xk < eps or xk > L - eps: continue - wx = dx/2 if (i == 0 or i == nx-1) else dx - Tx, Ty = traction(xk, yk, 0., +1., E, L) - F[k, 0] += Tx * wx; F[k, 1] += Ty * wx - return F # shape (N, 2) - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nx=10, ny=10, - with_visual=True, dim="2d"): - - tmpl = _dim_template(dim) - - rootNode.addObject("RequiredPlugin", pluginName="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - tris_np = _TriElement.get_triangles(nx, ny) - - Beam = rootNode.addChild("Beam") - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject('NewtonRaphsonSolver', - name="newtonSolver", - maxNbIterationsNewton=1, - absoluteResidualStoppingThreshold=1e-10, - printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), showObject=with_visual, - showObjectScale=0.005*L) - - Beam.addObject("TriangleSetTopologyContainer", name="topology", - triangles=tris_np.tolist()) - Beam.addObject("TriangleSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", name="FEM", template=tmpl, - youngModulus=E, poissonRatio=0.0, topology="@topology") - - e = 1e-4 * L / max(nx-1, ny-1) - if dim == "3d": - box_left = [-e, -e, -e, e, L+e, e] - box_right = [L-e, -e, -e, L+e, L+e, e] - else: - box_left = [-e, -e, 0, e, L+e, 0] - box_right = [L-e, -e, 0, L+e, L+e, 0] - - Beam.addObject("BoxROI", name="box_left", template=tmpl, - box=box_left, drawBoxes=with_visual) - Beam.addObject("FixedProjectiveConstraint", name="fix_left", template=tmpl, - indices="@box_left.indices") - Beam.addObject("BoxROI", name="box_right", template=tmpl, - box=box_right, drawBoxes=with_visual) - Beam.addObject("FixedProjectiveConstraint", name="fix_right", template=tmpl, - indices="@box_right.indices") - - # Build nodal force vector; pad with z=0 column for Vec3d - F_xy = _TriElement.compute_nodal_forces(nodes_2d, L, E, nx, ny) - if dim == "3d": - F_all = np.hstack([F_xy, np.zeros((len(F_xy), 1))]) - else: - F_all = F_xy - - idx = " ".join(str(k) for k in range(len(nodes_2d))) - if dim == "3d": - frc = " ".join(f"{F_all[k,0]} {F_all[k,1]} {F_all[k,2]}" - for k in range(len(nodes_2d))) - else: - frc = " ".join(f"{F_all[k,0]} {F_all[k,1]}" - for k in range(len(nodes_2d))) - - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - return dofs, nodes_2d, tris_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, tris_np): - """L2 error using centroid quadrature.""" - dx, dy = L / (nx - 1), L / (ny - 1) - area = dx * dy / 2.0 - xy = nodes_2d[:, :2] - err2 = 0.0 - for i0, i1, i2 in tris_np: - xc = (xy[i0,0] + xy[i1,0] + xy[i2,0]) / 3 - yc = (xy[i0,1] + xy[i1,1] + xy[i2,1]) / 3 - ux_c = (ux[i0] + ux[i1] + ux[i2]) / 3 - uy_c = (uy[i0] + uy[i1] + uy[i2]) / 3 - err2 += ((ux_c - ux_mms(xc,yc,L))**2 + - (uy_c - uy_mms(xc,yc,L))**2) * area - return np.sqrt(err2) - - @staticmethod - def _grad_p1(nodes_2d, u, tri): - xy = nodes_2d[:, :2] - i0, i1, i2 = tri - x0, y0 = xy[i0]; x1, y1 = xy[i1]; x2, y2 = xy[i2] - A2 = (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0) - dphi0_dx = (y1-y2)/A2; dphi0_dy = (x2-x1)/A2 - dphi1_dx = (y2-y0)/A2; dphi1_dy = (x0-x2)/A2 - dphi2_dx = (y0-y1)/A2; dphi2_dy = (x1-x0)/A2 - dudx = u[i0]*dphi0_dx + u[i1]*dphi1_dx + u[i2]*dphi2_dx - dudy = u[i0]*dphi0_dy + u[i1]*dphi1_dy + u[i2]*dphi2_dy - return dudx, dudy, abs(A2)/2.0 - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, tris_np, **_ignored): - err2 = 0.0 - for tri in tris_np: - i0, i1, i2 = tri - xy = nodes_2d[:, :2] - xc = (xy[i0,0] + xy[i1,0] + xy[i2,0]) / 3 - yc = (xy[i0,1] + xy[i1,1] + xy[i2,1]) / 3 - dux_dx_h, dux_dy_h, area = _TriElement._grad_p1(nodes_2d, ux, tri) - duy_dx_h, duy_dy_h, _ = _TriElement._grad_p1(nodes_2d, uy, tri) - err2 += ( - (dux_dx_h - dux_dx(xc,yc,L))**2 + - (dux_dy_h - dux_dy(xc,yc,L))**2 + - (duy_dx_h - duy_dx(xc,yc,L))**2 + - (duy_dy_h - duy_dy(xc,yc,L))**2 - ) * area - return np.sqrt(err2) - - -# ───────────────────────────────────────────────────────────────────────────── -# Simulation runner -# ───────────────────────────────────────────────────────────────────────────── - -def run_simulation(elem, L, E, nx, ny, dim="2d"): - root = Sofa.Core.Node("root") - dofs, nodes_2d, conn = elem.create_sofa_scene( - root, L=L, E=E, nx=nx, ny=ny, with_visual=False, dim=dim - ) - Sofa.Simulation.init(root) - pos0 = dofs.position.array().copy() - Sofa.Simulation.animate(root, root.dt.value) - pos1 = dofs.position.array().copy() - Sofa.Simulation.unload(root) - # Extract x and y displacement columns (works for both Vec2d and Vec3d) - ux = pos1[:, 0] - pos0[:, 0] - uy = pos1[:, 1] - pos0[:, 1] - return nodes_2d, conn, ux, uy - - -# ───────────────────────────────────────────────────────────────────────────── -# Convergence study -# ───────────────────────────────────────────────────────────────────────────── - -def convergence_study(elem_specs, L, E, nx_values, dim="2d", results_dir=None): - """ - Run convergence study for each element type in elem_specs. - elem_specs : list of dicts with keys 'elem', 'label', 'marker', 'color' - dim : "2d" or "3d" - """ - if results_dir is None: - results_dir = RESULTS_DIR - os.makedirs(results_dir, exist_ok=True) - - fig_l2, ax_l2 = plt.subplots(figsize=(7, 5)) - fig_h1, ax_h1 = plt.subplots(figsize=(7, 5)) - - for spec in elem_specs: - elem, label = spec["elem"], spec["label"] - marker, color = spec.get("marker", "o"), spec.get("color", None) - - hs, errs_l2, errs_h1 = [], [], [] - hdr = (f"{'nx':>5} | {'h':>8} | {'L2':>14} | {'ord_L2':>7} " - f"| {'H1':>14} | {'ord_H1':>7}") - tag = label.replace(" ", "_") - txt_path = os.path.join(results_dir, f"convergence_{tag}_{dim}.txt") - - print(f"\n── Convergence {label} [{dim.upper()}] ──\n{hdr}") - with open(txt_path, "w") as f: - f.write(f"Convergence MMS 2D — {label} [{dim.upper()}]\n{hdr}\n") - for k, nx in enumerate(nx_values): - ny = nx - h = L / (nx - 1) - nodes_2d, conn, ux, uy = run_simulation(elem, L, E, nx, ny, dim=dim) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - hs.append(h); errs_l2.append(l2); errs_h1.append(h1) - - ord_l2 = (f"{np.log(l2/errs_l2[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else "") - ord_h1 = (f"{np.log(h1/errs_h1[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else "") - line = (f"{nx:5d} | {h:8.4f} | {l2:14.6e} | {ord_l2:>7} " - f"| {h1:14.6e} | {ord_h1:>7}") - print(line); f.write(line + "\n") - - hs_a, l2_a, h1_a = np.array(hs), np.array(errs_l2), np.array(errs_h1) - kw = dict(lw=2, ms=7, color=color) - ax_l2.loglog(hs_a, l2_a, f"{marker}-", label=f"{label} [{dim}]", **kw) - ax_h1.loglog(hs_a, h1_a, f"{marker}--", label=f"{label} [{dim}]", **kw) - - # Reference slopes - h_ref = np.array([hs_a[0], hs_a[-1]]) - ax_l2.loglog(h_ref, l2_a[0]*(h_ref/hs_a[0])**2, ":", color="gray", - lw=1.2, label="O(h²)") - ax_h1.loglog(h_ref, h1_a[0]*(h_ref/hs_a[0])**1, ":", color="gray", - lw=1.2, label="O(h¹)") - - for ax, ylabel, title, fname in [ - (ax_l2, "L² error", f"L² convergence MMS 2D [{dim}]", - f"convergence_L2_{dim}.png"), - (ax_h1, "H¹ semi-norm", f"H¹ convergence MMS 2D [{dim}]", - f"convergence_H1_{dim}.png"), - ]: - ax.set_xlabel("h"); ax.set_ylabel(ylabel); ax.set_title(title) - ax.legend(); ax.grid(True, alpha=0.3, which="both") - - fig_l2.tight_layout(); fig_h1.tight_layout() - fig_l2.savefig(os.path.join(results_dir, f"convergence_L2_{dim}.png"), dpi=150) - fig_h1.savefig(os.path.join(results_dir, f"convergence_H1_{dim}.png"), dpi=150) - plt.close(fig_l2); plt.close(fig_h1) - - -# ───────────────────────────────────────────────────────────────────────────── -# Point simulation (solution plots) -# ───────────────────────────────────────────────────────────────────────────── - -def simulation_ponctuelle(elem, L, E, nx, ny, dim="2d", results_dir=None): - if results_dir is None: - results_dir = RESULTS_DIR - os.makedirs(results_dir, exist_ok=True) - - nodes_2d, conn, ux, uy = run_simulation(elem, L, E, nx, ny, dim=dim) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - - xy = nodes_2d[:, :2] - ux_ref = ux_mms(xy[:, 0], xy[:, 1], L) - uy_ref = uy_mms(xy[:, 0], xy[:, 1], L) - - label = elem.LABEL - print(f"\n[{label}] dim={dim} nx={nx} ny={ny} L={L}") - print(f" L2 = {l2:.4e}") - print(f" H1 semi-norm = {h1:.4e}") - print(f" max|ux-ux_mms| = {np.max(np.abs(ux - ux_ref)):.4e}") - print(f" max|uy-uy_mms| = {np.max(np.abs(uy - uy_ref)):.4e}") - - # 1-D profile along mid-height line - mid_j = (ny - 1) // 2 - sl = slice(mid_j * nx, mid_j * nx + nx) - yc = xy[mid_j * nx, 1] - xf = np.linspace(0, L, 300) - - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) - for ax, u_sofa, u_fn, lbl, fmt in zip( - axes, - [ux[sl], uy[sl]], - [lambda x: ux_mms(x, yc, L), lambda x: uy_mms(x, yc, L)], - [r"$u_x$", r"$u_y$"], - ["o-", "s-"], - ): - ax.plot(xy[sl, 0], u_sofa, fmt, color="tab:green", - label=f"SOFA {label} [{dim}]", ms=5) - ax.plot(xf, u_fn(xf), "--", color="tab:blue", label="MMS exact") - ax.set_xlabel("x"); ax.set_ylabel(lbl) - ax.legend(); ax.grid(True, alpha=0.3) - plt.suptitle(f"MMS 2D — {label} [{dim}] nx={nx} |L2={l2:.2e} H1={h1:.2e}") - plt.tight_layout() - tag = label.replace(" ", "_") - plt.savefig(os.path.join(results_dir, - f"solution_{tag}_{dim}_nx{nx}.png"), dpi=150) - plt.close() - - # 2-D colour maps - tris_plot = elem.to_triangles(conn) - x, y = xy[:, 0], xy[:, 1] - fig, axes = plt.subplots(2, 3, figsize=(15, 8)) - for ax, data, title, cmap in [ - (axes[0,0], ux, r"$u_x$ SOFA", "gray"), - (axes[0,1], ux_ref, r"$u_x$ MMS", "gray"), - (axes[0,2], np.abs(ux-ux_ref), r"$|u_x - u_x^{MMS}|$", "gray_r"), - (axes[1,0], uy, r"$u_y$ SOFA", "gray"), - (axes[1,1], uy_ref, r"$u_y$ MMS", "gray"), - (axes[1,2], np.abs(uy-uy_ref), r"$|u_y - u_y^{MMS}|$", "gray_r"), - ]: - tc = ax.tricontourf(x, y, tris_plot.tolist(), data, levels=20, cmap=cmap) - ax.triplot(x, y, tris_plot.tolist(), "k-", lw=0.3, alpha=0.4) - plt.colorbar(tc, ax=ax, shrink=0.8) - ax.set_title(title); ax.set_aspect("equal") - ax.set_xlabel("x"); ax.set_ylabel("y") - plt.suptitle(f"Champs 2D — {label} [{dim}] nx={nx}") - plt.tight_layout() - plt.savefig(os.path.join(results_dir, - f"champs2D_{tag}_{dim}_nx{nx}.png"), dpi=150) - plt.close() - - -# ───────────────────────────────────────────────────────────────────────────── -# Element instances -# ───────────────────────────────────────────────────────────────────────────── - -element_quad = _QuadElement() -element_tri = _TriElement() - - -if __name__ == "__main__": - - - # DIM = "2d" Vec2d template (plan) - # DIM = "3d" Vec3d template (z=0) - - DIM = "2d" - - L, E = 1.0, 1e6 - - nx_vals = [10,20,30,40,50,60] - - - - specs = [ - {"elem": element_quad, "label": "Q1 quad", "marker": "o", "color": "C0"}, - {"elem": element_tri, "label": "P1 tri", "marker": "s", "color": "C1"}, - ] - - convergence_study(specs, L, E, nx_vals, dim=DIM) - simulation_ponctuelle(element_quad, L, E, nx=20, ny=20, dim=DIM) \ No newline at end of file diff --git a/examples/Freefem/MMS/2D/params.json b/examples/Freefem/MMS/2D/params.json new file mode 100644 index 00000000..66b001d2 --- /dev/null +++ b/examples/Freefem/MMS/2D/params.json @@ -0,0 +1,17 @@ +{ + "length": 1.0, + "youngModulus": 1e6, + "reference": { + "nx": 20, + "nu": 0.0, + "dim": "2d" + }, + "convergence": { + "nu_values": [0.0, 0.3, 0.49], + "dim_values": ["2d", "3d"], + "nx_values": { + "cubic": [10, 20, 30, 40, 50], + "trigonometric": [10, 20, 30, 40, 50] + } + } +} diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py new file mode 100644 index 00000000..724b73d3 --- /dev/null +++ b/examples/Freefem/MMS/2D/run_convergence.py @@ -0,0 +1,163 @@ +"""Run the mesh-refinement convergence study for every 2D MMS case. + +Loops cases × elements × dim × nu × nx and writes per-(case, element, dim, nu) +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 beam import ( + RESULTS_DIR, + load_params, + solve_beam, + element_quad, + element_tri, +) +from plot_utils import annotate_convergence_rates + + +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 + 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' + dim : "2d" (plane stress) or "3d" (plane strain) + """ + 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: + elem, label = spec["elem"], spec["label"] + 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) + + rate_l2 = (f"{np.log(l2 / l2s[-1]) / np.log(h / hs[-1]):.2f}" + if k > 0 else "") + rate_h1 = (f"{np.log(h1 / h1s[-1]) / np.log(h / hs[-1]):.2f}" + if k > 0 else "") + print(f"{nx:5d} | {h:10.4f} | {l2:14.6e} | {rate_l2:>7} " + f"| {h1:14.6e} | {rate_h1:>7}", flush=True) + rows.append({"nx": nx, "h": h, + "L2": l2, "rate_L2": rate_l2, + "H1": h1, "rate_H1": rate_h1}) + hs.append(h); l2s.append(l2); h1s.append(h1) + + write_convergence_table(stem, rows) + plot_series.append({"label": f"{label} L²", + "errors": l2s, "style": spec["l2_style"]}) + plot_series.append({"label": f"{label} H¹", + "errors": h1s, "style": spec["h1_style"]}) + hs_ref = hs + + first_label = elem_specs[0]["label"] + ref_lines = [ + {"slope": 2, "anchor": f"{first_label} L²", "label": "O(h²)"}, + {"slope": 1, "anchor": f"{first_label} H¹", "label": "O(h¹)"}, + ] + title = f"Convergence — {mms.name} [{dim} / {hyp}] nu={nu}" + plot_convergence(f"convergence_{mms.name}_{dim}_nu{nu}", + hs_ref, plot_series, title=title, ref_lines=ref_lines) + + +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", ref_lines=None): + """ + Save log-log convergence plot to results/.png. + + series : list of {"label", "errors", "style"?} dicts + ref_lines : optional list of {"slope", "anchor", "label"?} dicts; + "anchor" names the series whose first point anchors the line. + 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)) + anchors = {} + 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) + anchors[s["label"]] = e_arr + for rl in (ref_lines or []): + anchor = anchors.get(rl["anchor"]) + if anchor is None: + continue + slope = rl["slope"] + label = rl.get("label", f"O(h^{slope})") + h_ref = np.array([h_arr[0], h_arr[-1]]) + ax.loglog(h_ref, anchor[0] * (h_ref / h_arr[0])**slope, ":", + color="gray", lw=1.2, label=label) + ax.set_xlabel("h") + ax.set_ylabel(ylabel) + ax.set_title(title) + ax.legend() + ax.grid(True, alpha=0.3, which="both") + fig.tight_layout() + fig.savefig(os.path.join(RESULTS_DIR, f"{stem}.png"), dpi=150) + plt.close(fig) + + +if __name__ == "__main__": + cfg = load_params() + L = cfg["length"] + E = cfg["youngModulus"] + conv = cfg["convergence"] + + specs = [ + {"elem": element_quad, "label": "Q1 quad", + "l2_style": "bo-", "h1_style": "rs--"}, + {"elem": element_tri, "label": "P1 tri", + "l2_style": "b^-", "h1_style": "rD--"}, + ] + + # dim="2d" → Vec2d / plane stress; dim="3d" → Vec3d / plane strain + for mms in (cubic_mms, trig_mms): + nx_vals = conv["nx_values"][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/2D/trigonometric-cubic-mms.py b/examples/Freefem/MMS/2D/trigonometric-cubic-mms.py deleted file mode 100644 index 6fd81824..00000000 --- a/examples/Freefem/MMS/2D/trigonometric-cubic-mms.py +++ /dev/null @@ -1,836 +0,0 @@ -""" -MMS 2D - elasticity linear -2 MMS in this code : -cubic : ux = x(L-x)(x+y), uy = y(L-y)(y-x) -trigo : ux = sin(πx/L) cos(πy/L), uy = cos(πx/L) sin(πy/L) - -Q1 & P1 elements -nu = Poisson ratio -dim="2d" ==> Vec2d ==> plane stress \lambda = E*nu/(1-nu^2) -dim="3d" ==> Vec3d ==> plane strain \lambda = E* nu/((1+nu)(1-2*nu)) -""" - -import numpy as np -import matplotlib.pyplot as plt -import os - -import Sofa -import Sofa.Core -import Sofa.Simulation -import SofaRuntime - -# =========== Choose the MMS ============================ -MMS_TYPE = "cubic" # "cubic" ou "trigonometric" - -RESULTS_DIR = f"results_{MMS_TYPE}_q1p1" -os.makedirs(RESULTS_DIR, exist_ok=True) - -GAUSS_PTS = np.array([-1.0 / np.sqrt(3), 1.0 / np.sqrt(3)]) -GAUSS_WTS = np.array([1.0, 1.0]) - -# ============== LAME ==================== -def lame(E, nu, dim="2d"): - """ - dim="2d" ===> plane stress : \lambda = E*nu/(1-nu^²) - dim="3d" ==> Vec3d ==> plane strain \lambda = E* nu/((1+nu)(1-2*nu)) - - """ - if dim == "2d": - lam = E * nu / (1 - nu**2) - else: - lam = E * nu / ((1 + nu) * (1 - 2 * nu)) - mu = E / (2 * (1 + nu)) - return lam, mu - -# ==================== Cubic MMS ================================= -def ux_mms_cubic(x, y, L): - return x * (L - x) * (x + y) - -def uy_mms_cubic(x, y, L): - return y * (L - y) * (y - x) - -# --------------- Derivatives cubic --------------- -def dux_dx_cubic(x, y, L): - return (L - 2*x) * (x + y) + x * (L - x) - -def dux_dy_cubic(x, y, L): - return x * (L - x) - -def duy_dx_cubic(x, y, L): - return -y * (L - y) - -def duy_dy_cubic(x, y, L): - return (L - 2*y) * (y - x) + y * (L - y) - -def d2ux_dxx_cubic(x, y, L): - return 2*L - 6*x - 2*y - -def d2ux_dxy_cubic(x, y, L): - return L - 2*x - -def d2ux_dyy_cubic(x, y, L): - return np.zeros_like(np.asarray(x, float)) - -def d2uy_dxx_cubic(x, y, L): - return np.zeros_like(np.asarray(x, float)) - -def d2uy_dxy_cubic(x, y, L): - return -(L - 2*y) - -def d2uy_dyy_cubic(x, y, L): - return 2*L - 6*y + 2*x - -# --------------- Body force cubic --------------- -def fx_body_cubic(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return -((lam + 2*mu) * d2ux_dxx_cubic(x, y, L) - + lam * d2uy_dxy_cubic(x, y, L) - + mu * (d2ux_dyy_cubic(x, y, L) + d2uy_dxy_cubic(x, y, L))) - -def fy_body_cubic(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return -(mu * (d2ux_dxy_cubic(x, y, L) + d2uy_dxx_cubic(x, y, L)) - + lam * d2ux_dxy_cubic(x, y, L) - + (lam + 2*mu) * d2uy_dyy_cubic(x, y, L)) - -# --------------- Stress cubic --------------- -def sigma_xx_cubic(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return (lam + 2*mu) * dux_dx_cubic(x, y, L) + lam * duy_dy_cubic(x, y, L) - -def sigma_yy_cubic(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return lam * dux_dx_cubic(x, y, L) + (lam + 2*mu) * duy_dy_cubic(x, y, L) - -def sigma_xy_cubic(x, y, E, nu, L, dim="2d"): - _, mu = lame(E, nu, dim) - return mu * (dux_dy_cubic(x, y, L) + duy_dx_cubic(x, y, L)) - -def traction_cubic(x, y, nx_c, ny_c, E, nu, L, dim="2d"): - sxx = sigma_xx_cubic(x, y, E, nu, L, dim) - syy = sigma_yy_cubic(x, y, E, nu, L, dim) - sxy = sigma_xy_cubic(x, y, E, nu, L, dim) - return sxx*nx_c + sxy*ny_c, sxy*nx_c + syy*ny_c - -# ==================== Trigonometric MMS =================================== -def ux_mms_trigo(x, y, L): - return np.sin(np.pi * x / L) * np.cos(np.pi * y / L) - -def uy_mms_trigo(x, y, L): - return np.cos(np.pi * x / L) * np.sin(np.pi * y / L) - -# --------------- Trigo derivatives --------------- -def dux_dx_trigo(x, y, L): - return (np.pi / L) * np.cos(np.pi * x / L) * np.cos(np.pi * y / L) - -def dux_dy_trigo(x, y, L): - return -(np.pi / L) * np.sin(np.pi * x / L) * np.sin(np.pi * y / L) - -def duy_dx_trigo(x, y, L): - return -(np.pi / L) * np.sin(np.pi * x / L) * np.sin(np.pi * y / L) - -def duy_dy_trigo(x, y, L): - return (np.pi / L) * np.cos(np.pi * x / L) * np.cos(np.pi * y / L) - -def d2ux_dxx_trigo(x, y, L): - return -(np.pi / L)**2 * ux_mms_trigo(x, y, L) - -def d2ux_dyy_trigo(x, y, L): - return -(np.pi / L)**2 * ux_mms_trigo(x, y, L) - -def d2ux_dxy_trigo(x, y, L): - return -(np.pi / L)**2 * np.cos(np.pi * x / L) * np.sin(np.pi * y / L) - -def d2uy_dxx_trigo(x, y, L): - return -(np.pi / L)**2 * uy_mms_trigo(x, y, L) - -def d2uy_dyy_trigo(x, y, L): - return -(np.pi / L)**2 * uy_mms_trigo(x, y, L) - -def d2uy_dxy_trigo(x, y, L): - return -(np.pi / L)**2 * np.sin(np.pi * x / L) * np.cos(np.pi * y / L) - -# --------------- Trigo body force --------------- -def fx_body_trigo(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - ddiv_dx = d2ux_dxx_trigo(x, y, L) + d2uy_dxy_trigo(x, y, L) - laplace_ux = d2ux_dxx_trigo(x, y, L) + d2ux_dyy_trigo(x, y, L) - return -(lam + mu) * ddiv_dx - mu * laplace_ux - -def fy_body_trigo(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - ddiv_dy = d2ux_dxy_trigo(x, y, L) + d2uy_dyy_trigo(x, y, L) - laplace_uy = d2uy_dxx_trigo(x, y, L) + d2uy_dyy_trigo(x, y, L) - return -(lam + mu) * ddiv_dy - mu * laplace_uy - -# --------------- Stress trigo --------------- -def sigma_xx_trigo(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return (lam + 2*mu) * dux_dx_trigo(x, y, L) + lam * duy_dy_trigo(x, y, L) - -def sigma_yy_trigo(x, y, E, nu, L, dim="2d"): - lam, mu = lame(E, nu, dim) - return lam * dux_dx_trigo(x, y, L) + (lam + 2*mu) * duy_dy_trigo(x, y, L) - -def sigma_xy_trigo(x, y, E, nu, L, dim="2d"): - _, mu = lame(E, nu, dim) - return mu * (dux_dy_trigo(x, y, L) + duy_dx_trigo(x, y, L)) - -def traction_trigo(x, y, nx_c, ny_c, E, nu, L, dim="2d"): - sxx = sigma_xx_trigo(x, y, E, nu, L, dim) - syy = sigma_yy_trigo(x, y, E, nu, L, dim) - sxy = sigma_xy_trigo(x, y, E, nu, L, dim) - return sxx*nx_c + sxy*ny_c, sxy*nx_c + syy*ny_c - -# ==================== MMS Selection ==================== -if MMS_TYPE == "cubic": - ux_mms = ux_mms_cubic - uy_mms = uy_mms_cubic - dux_dx = dux_dx_cubic - dux_dy = dux_dy_cubic - duy_dx = duy_dx_cubic - duy_dy = duy_dy_cubic - fx_body = fx_body_cubic - fy_body = fy_body_cubic - sigma_xx = sigma_xx_cubic - sigma_yy = sigma_yy_cubic - sigma_xy = sigma_xy_cubic - traction = traction_cubic - print(" cubic MMS activated") -elif MMS_TYPE == "trigonometric": - ux_mms = ux_mms_trigo - uy_mms = uy_mms_trigo - dux_dx = dux_dx_trigo - dux_dy = dux_dy_trigo - duy_dx = duy_dx_trigo - duy_dy = duy_dy_trigo - fx_body = fx_body_trigo - fy_body = fy_body_trigo - sigma_xx = sigma_xx_trigo - sigma_yy = sigma_yy_trigo - sigma_xy = sigma_xy_trigo - traction = traction_trigo - print(" trigonometric MMS activated") -else: - raise ValueError(f"MMS_TYPE inconnu: {MMS_TYPE}") - -# ==================== MESH ==================== -def get_nodes_2d(L, nx, ny, dim="2d"): - dx, dy = L / (nx - 1), L / (ny - 1) - if dim == "3d": - pts = [[i*dx, j*dy, 0.0] for j in range(ny) for i in range(nx)] - else: - pts = [[i*dx, j*dy] for j in range(ny) for i in range(nx)] - return np.array(pts, dtype=float) - -def _sofa_template(dim): - return "Vec3d" if dim == "3d" else "Vec2d" - -# ==================== DIRICHLET ==================== -def _add_dirichlet(Beam, nodes_2d, L, dim, with_visual): - tmpl = _sofa_template(dim) - xy = nodes_2d[:, :2] - eps = 1e-10 - - idx_ux0 = [k for k, (xk, _) in enumerate(xy) if xk < eps or xk > L - eps] - idx_uy0 = [k for k, (_, yk) in enumerate(xy) if yk < eps or yk > L - eps] - - if dim == "2d": - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_ux", template=tmpl, - indices=" ".join(map(str, idx_ux0)), - fixedDirections="1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_uy", template=tmpl, - indices=" ".join(map(str, idx_uy0)), - fixedDirections="0 1") - else: - all_idx = " ".join(str(k) for k in range(len(nodes_2d))) - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_ux", template=tmpl, - indices=" ".join(map(str, idx_ux0)), - fixedDirections="1 0 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_uy", template=tmpl, - indices=" ".join(map(str, idx_uy0)), - fixedDirections="0 1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_z", template=tmpl, - indices=all_idx, - fixedDirections="0 0 1") - -def _pack_forces(F_xy, dim, n_nodes): - idx = " ".join(str(k) for k in range(n_nodes)) - if dim == "3d": - frc = " ".join(f"{F_xy[k,0]} {F_xy[k,1]} 0.0" for k in range(n_nodes)) - else: - frc = " ".join(f"{F_xy[k,0]} {F_xy[k,1]}" for k in range(n_nodes)) - return idx, frc - -# ==================== Q1 element ==================== -class _QuadElement: - LABEL = "Q1 quad" - - @staticmethod - def get_connectivity(nx, ny): - quads = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - quads.append([k00, k10, k11, k01]) - return np.array(quads) - - @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) - - @staticmethod - def _shape(xi, eta): - 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)]) - return N, dN_dxi, dN_deta - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim="2d"): - quads = _QuadElement.get_connectivity(nx, ny) - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - - - for quad in quads: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(xi, eta) - J = np.array([[dN_dxi @ xe, dN_dxi @ ye], - [dN_deta @ xe, dN_deta @ ye]]) - detJ = np.linalg.det(J) - xg, yg = N @ xe, N @ ye - - fx = fx_body(xg, yg, E, nu, L, dim) - fy = fy_body(xg, yg, E, nu, L, dim) - w = wi * wj * detJ - for a, node in enumerate(quad): - F[node, 0] += N[a] * fx * w - F[node, 1] += N[a] * fy * w - - # Neumann x=0 - for j in range(ny - 1): - k0, k1 = j * nx, (j + 1) * nx - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5 * (xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, -1., 0., E, nu, L, dim) - N0, N1 = 0.5*(1 - xi), 0.5*(1 + xi) - w = wi * Le / 2.0 - F[k0, 0] += N0 * Tx * w; F[k0, 1] += N0 * Ty * w - F[k1, 0] += N1 * Tx * w; F[k1, 1] += N1 * Ty * w - - # Neumann x=L - for j in range(ny - 1): - k0 = j * nx + (nx - 1) - k1 = (j + 1) * nx + (nx - 1) - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5 * (xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, +1., 0., E, nu, L, dim) - N0, N1 = 0.5*(1 - xi), 0.5*(1 + xi) - w = wi * Le / 2.0 - F[k0, 0] += N0 * Tx * w; F[k0, 1] += N0 * Ty * w - F[k1, 0] += N1 * Tx * w; F[k1, 1] += N1 * Ty * w - - # Neumann y=0 - for i in range(nx - 1): - k0, k1 = i, i + 1 - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5 * (xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, 0., -1., E, nu, L, dim) - N0, N1 = 0.5*(1 - xi), 0.5*(1 + xi) - w = wi * Le / 2.0 - F[k0, 0] += N0 * Tx * w; F[k0, 1] += N0 * Ty * w - F[k1, 0] += N1 * Tx * w; F[k1, 1] += N1 * Ty * w - - # Neumann y=L - for i in range(nx - 1): - k0 = (ny - 1) * nx + i - k1 = (ny - 1) * nx + i + 1 - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5 * (xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction(xg, yg, 0., +1., E, nu, L, dim) - N0, N1 = 0.5*(1 - xi), 0.5*(1 + xi) - w = wi * Le / 2.0 - F[k0, 0] += N0 * Tx * w; F[k0, 1] += N0 * Ty * w - F[k1, 0] += N1 * Tx * w; F[k1, 1] += N1 * Ty * w - - return F - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nu=0.3, - nx=10, ny=10, with_visual=True, dim="2d"): - tmpl = _sofa_template(dim) - - rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - quads_np = _QuadElement.get_connectivity(nx, ny) - - Beam = rootNode.addChild("Beam") - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), - showObject=with_visual, showObjectScale=0.005*L) - - Beam.addObject("QuadSetTopologyContainer", - name="topology", quads=quads_np.tolist()) - Beam.addObject("QuadSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", - name="FEM", template=tmpl, - youngModulus=E, poissonRatio=nu, topology="@topology") - - _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) - - - F_xy = _QuadElement.compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim) - idx, frc = _pack_forces(F_xy, dim, len(nodes_2d)) - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - - return dofs, nodes_2d, quads_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, conn): - xy, err2 = nodes_2d[:, :2], 0.0 - for quad in conn: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(xi, eta) - J = np.array([[dN_dxi @ xe, dN_dxi @ ye], - [dN_deta @ xe, dN_deta @ ye]]) - detJ = np.linalg.det(J) - xg, yg = N @ xe, N @ ye - ux_h = N @ ux[quad] - uy_h = N @ uy[quad] - err2 += ((ux_h - ux_mms(xg, yg, L))**2 - + (uy_h - uy_mms(xg, yg, L))**2) * wi * wj * detJ - return np.sqrt(err2) - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, conn, **_): - xy, err2 = nodes_2d[:, :2], 0.0 - for quad in conn: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(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 - err2 += ( - (dN_dx @ ux[quad] - dux_dx(xg, yg, L))**2 + - (dN_dy @ ux[quad] - dux_dy(xg, yg, L))**2 + - (dN_dx @ uy[quad] - duy_dx(xg, yg, L))**2 + - (dN_dy @ uy[quad] - duy_dy(xg, yg, L))**2 - ) * wi * wj * detJ - return np.sqrt(err2) - -# ==================== P1 elements ==================== -class _TriElement: - LABEL = "P1 tri" - - @staticmethod - def get_connectivity(nx, ny): - tris = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - tris += [[k00, k10, k11], [k00, k11, k01]] - return np.array(tris) - - @staticmethod - def to_triangles(conn): - return conn - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim="2d"): - dx, dy = L / (nx - 1), L / (ny - 1) - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - eps = 1e-10 - - - for k, (xk, yk) in enumerate(xy): - i = round(xk / dx) - j = round(yk / dy) - wx = dx/2 if (i == 0 or i == nx-1) else dx - wy = dy/2 if (j == 0 or j == ny-1) else dy - # dim transmis - F[k, 0] += fx_body(xk, yk, E, nu, L, dim) * wx * wy - F[k, 1] += fy_body(xk, yk, E, nu, L, dim) * wx * wy - - # Neumann x=0 - for j in range(ny): - k = j * nx - xk, yk = xy[k] - if yk < eps or yk > L - eps: - continue - wy = dy/2 if (j == 0 or j == ny-1) else dy - Tx, Ty = traction(xk, yk, -1., 0., E, nu, L, dim) - F[k, 0] += Tx * wy; F[k, 1] += Ty * wy - - # Neumann x=L - for j in range(ny): - k = j * nx + (nx - 1) - xk, yk = xy[k] - if yk < eps or yk > L - eps: - continue - wy = dy/2 if (j == 0 or j == ny-1) else dy - Tx, Ty = traction(xk, yk, +1., 0., E, nu, L, dim) - F[k, 0] += Tx * wy; F[k, 1] += Ty * wy - - # Neumann y=0 - for i in range(nx): - k = i - xk, yk = xy[k] - if xk < eps or xk > L - eps: - continue - wx = dx/2 if (i == 0 or i == nx-1) else dx - Tx, Ty = traction(xk, yk, 0., -1., E, nu, L, dim) - F[k, 0] += Tx * wx; F[k, 1] += Ty * wx - - # Neumann y=L - for i in range(nx): - k = (ny - 1) * nx + i - xk, yk = xy[k] - if xk < eps or xk > L - eps: - continue - wx = dx/2 if (i == 0 or i == nx-1) else dx - Tx, Ty = traction(xk, yk, 0., +1., E, nu, L, dim) - F[k, 0] += Tx * wx; F[k, 1] += Ty * wx - - return F - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nu=0.3, - nx=10, ny=10, with_visual=True, dim="2d"): - tmpl = _sofa_template(dim) - - rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - tris_np = _TriElement.get_connectivity(nx, ny) - - Beam = rootNode.addChild("Beam") - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), - showObject=with_visual, showObjectScale=0.005*L) - - Beam.addObject("TriangleSetTopologyContainer", - name="topology", triangles=tris_np.tolist()) - Beam.addObject("TriangleSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", - name="FEM", template=tmpl, - youngModulus=E, poissonRatio=nu, topology="@topology") - - _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) - - - F_xy = _TriElement.compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim) - idx, frc = _pack_forces(F_xy, dim, len(nodes_2d)) - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - - return dofs, nodes_2d, tris_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, conn): - dx, dy = L / (nx - 1), L / (ny - 1) - area = dx * dy / 2.0 - xy = nodes_2d[:, :2] - err2 = 0.0 - for i0, i1, i2 in conn: - xc = (xy[i0,0] + xy[i1,0] + xy[i2,0]) / 3 - yc = (xy[i0,1] + xy[i1,1] + xy[i2,1]) / 3 - ux_c = (ux[i0] + ux[i1] + ux[i2]) / 3 - uy_c = (uy[i0] + uy[i1] + uy[i2]) / 3 - err2 += ((ux_c - ux_mms(xc, yc, L))**2 - + (uy_c - uy_mms(xc, yc, L))**2) * area - return np.sqrt(err2) - - @staticmethod - def _grad_p1(xy, u, tri): - i0, i1, i2 = tri - x0, y0 = xy[i0] - x1, y1 = xy[i1] - x2, y2 = xy[i2] - A2 = (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0) - dudx = (u[i0]*(y1-y2) + u[i1]*(y2-y0) + u[i2]*(y0-y1)) / A2 - dudy = (u[i0]*(x2-x1) + u[i1]*(x0-x2) + u[i2]*(x1-x0)) / A2 - return dudx, dudy, abs(A2) / 2.0 - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, conn, **_): - xy, err2 = nodes_2d[:, :2], 0.0 - for tri in conn: - i0, i1, i2 = tri - xc = (xy[i0,0] + xy[i1,0] + xy[i2,0]) / 3 - yc = (xy[i0,1] + xy[i1,1] + xy[i2,1]) / 3 - dux_dx_h, dux_dy_h, area = _TriElement._grad_p1(xy, ux, tri) - duy_dx_h, duy_dy_h, _ = _TriElement._grad_p1(xy, uy, tri) - err2 += ( - (dux_dx_h - dux_dx(xc, yc, L))**2 + - (dux_dy_h - dux_dy(xc, yc, L))**2 + - (duy_dx_h - duy_dx(xc, yc, L))**2 + - (duy_dy_h - duy_dy(xc, yc, L))**2 - ) * area - return np.sqrt(err2) - -# ==================== SIMULATION ==================== -def run_simulation(elem, L, E, nu, nx, ny, dim="2d"): - root = Sofa.Core.Node("root") - dofs, nodes_2d, conn = elem.create_sofa_scene( - root, L=L, E=E, nu=nu, nx=nx, ny=ny, with_visual=False, dim=dim - ) - Sofa.Simulation.init(root) - pos0 = dofs.position.array().copy() - Sofa.Simulation.animate(root, root.dt.value) - pos1 = dofs.position.array().copy() - Sofa.Simulation.unload(root) - ux = pos1[:, 0] - pos0[:, 0] - uy = pos1[:, 1] - pos0[:, 1] - return nodes_2d, conn, ux, uy - -# ==================== POST-TRAITEMENT ==================== -def simulation_ponctuelle(elem, L, E, nu, nx, ny, dim="2d", results_dir=RESULTS_DIR): - os.makedirs(results_dir, exist_ok=True) - nodes_2d, conn, ux, uy = run_simulation(elem, L, E, nu, nx, ny, dim=dim) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - - xy = nodes_2d[:, :2] - ux_ref = ux_mms(xy[:, 0], xy[:, 1], L) - uy_ref = uy_mms(xy[:, 0], xy[:, 1], L) - - hyp = "plane strain" if dim == "3d" else "plane stress" - print(f"\n[{elem.LABEL}] {MMS_TYPE.upper()} dim={dim} ({hyp}) nu={nu} nx={nx} L={L}") - print(f" L2 = {l2:.4e}") - print(f" H1 semi-norm = {h1:.4e}") - print(f" max|ux-ux_mms| = {np.max(np.abs(ux - ux_ref)):.4e}") - print(f" max|uy-uy_mms| = {np.max(np.abs(uy - uy_ref)):.4e}") - - mid_j = (ny - 1) // 2 - sl = slice(mid_j * nx, mid_j * nx + nx) - yc = xy[mid_j * nx, 1] - xf = np.linspace(0, L, 300) - - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) - for ax, u_sofa, u_fn, lbl, fmt in zip( - axes, - [ux[sl], uy[sl]], - [lambda x: ux_mms(x, yc, L), lambda x: uy_mms(x, yc, L)], - [r"$u_x$", r"$u_y$"], ["o-", "s-"], - ): - ax.plot(xy[sl, 0], u_sofa, fmt, color="tab:orange", - label=f"SOFA {elem.LABEL} [{dim}]", ms=5) - ax.plot(xf, u_fn(xf), "--", color="tab:green", label="MMS exact") - ax.set_xlabel("x") - ax.set_ylabel(lbl) - ax.legend() - ax.grid(True, alpha=0.3) - plt.suptitle( - f"MMS {MMS_TYPE.upper()} — {elem.LABEL} [{dim} / {hyp}] nu={nu} nx={nx}" - f" |L2={l2:.2e} H1={h1:.2e}" - ) - plt.tight_layout() - tag = elem.LABEL.replace(" ", "_") - plt.savefig(os.path.join(results_dir, - f"solution_{tag}_{dim}_nu{nu}_nx{nx}.png"), dpi=150) - plt.close() - - tris_plot = elem.to_triangles(conn) - x, y = xy[:, 0], xy[:, 1] - fig, axes = plt.subplots(2, 3, figsize=(15, 8)) - for ax, data, title, cmap in [ - (axes[0,0], ux, r"$u_x$ SOFA", "gray"), - (axes[0,1], ux_ref, r"$u_x$ MMS", "gray"), - (axes[0,2], np.abs(ux-ux_ref), r"$|u_x - u_x^{MMS}|$", "gray_r"), - (axes[1,0], uy, r"$u_y$ SOFA", "gray"), - (axes[1,1], uy_ref, r"$u_y$ MMS", "gray"), - (axes[1,2], np.abs(uy-uy_ref), r"$|u_y - u_y^{MMS}|$", "gray_r"), - ]: - tc = ax.tricontourf(x, y, tris_plot.tolist(), data, levels=20, cmap=cmap) - ax.triplot(x, y, tris_plot.tolist(), "k-", lw=0.3, alpha=0.4) - plt.colorbar(tc, ax=ax, shrink=0.8) - ax.set_title(title) - ax.set_aspect("equal") - ax.set_xlabel("x"); ax.set_ylabel("y") - plt.suptitle( - f"Champs 2D — {elem.LABEL} [{dim} / {hyp}] {MMS_TYPE.upper()} nu={nu} nx={nx}" - ) - plt.tight_layout() - plt.savefig(os.path.join(results_dir, - f"champs2D_{tag}_{dim}_nu{nu}_nx{nx}.png"), dpi=150) - plt.close() - -# ==================== CONVERGENCE study ==================== -def convergence_study(elem_specs, L, E, nu, nx_values, dim="2d", results_dir=RESULTS_DIR): - os.makedirs(results_dir, exist_ok=True) - hyp = "plane strain" if dim == "3d" else "plane stress" - fig_l2, ax_l2 = plt.subplots(figsize=(7, 5)) - fig_h1, ax_h1 = plt.subplots(figsize=(7, 5)) - - hs_last = l2_last = h1_last = None - - for spec in elem_specs: - elem = spec["elem"] - label = spec["label"] - marker = spec.get("marker", "o") - color = spec.get("color", None) - - hs, errs_l2, errs_h1 = [], [], [] - hdr = (f"{'nx':>5} | {'h':>8} | {'L2':>14} | {'ord_L2':>7}" - f" | {'H1':>14} | {'ord_H1':>7}") - tag = label.replace(" ", "_") - txt_path = os.path.join(results_dir, - f"convergence_{tag}_{dim}_nu{nu}.txt") - - print(f"\n── Convergence {label} [{dim} / {hyp}] {MMS_TYPE.upper()} nu={nu} ──\n{hdr}") - with open(txt_path, "w", encoding="utf-8") as f: - f.write(f"Convergence MMS {MMS_TYPE.upper()} — {label} [{dim} / {hyp}] nu={nu}\n{hdr}\n") - for k, nx in enumerate(nx_values): - ny = nx - h = L / (nx - 1) - nodes_2d, conn, ux, uy = run_simulation( - elem, L, E, nu, nx, ny, dim=dim - ) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - hs.append(h); errs_l2.append(l2); errs_h1.append(h1) - - ord_l2 = (f"{np.log(l2/errs_l2[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else " — ") - ord_h1 = (f"{np.log(h1/errs_h1[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else " — ") - line = (f"{nx:5d} | {h:8.4f} | {l2:14.6e} | {ord_l2:>7}" - f" | {h1:14.6e} | {ord_h1:>7}") - print(line) - f.write(line + "\n") - - hs_a = np.array(hs) - l2_a = np.array(errs_l2) - h1_a = np.array(errs_h1) - kw = dict(lw=2, ms=7, color=color) - ax_l2.loglog(hs_a, l2_a, f"{marker}-", - label=f"{label} [{dim}] {MMS_TYPE.upper()} nu={nu}", **kw) - ax_h1.loglog(hs_a, h1_a, f"{marker}--", - label=f"{label} [{dim}] {MMS_TYPE.upper()} nu={nu}", **kw) - hs_last, l2_last, h1_last = hs_a, l2_a, h1_a - - if hs_last is not None: - h_ref = np.array([hs_last[0], hs_last[-1]]) - ax_l2.loglog(h_ref, l2_last[0]*(h_ref/hs_last[0])**2, - ":", color="gray", lw=1.2, label="O(h²)") - ax_h1.loglog(h_ref, h1_last[0]*(h_ref/hs_last[0])**1, - ":", color="gray", lw=1.2, label="O(h¹)") - - for ax, ylabel, title in [ - (ax_l2, "L2 error", f"Convergence L2 — MMS {MMS_TYPE.upper()} [{dim} / {hyp}]"), - (ax_h1, "H1 semi-norm", f"Convergence H1 — MMS {MMS_TYPE.upper()} [{dim} / {hyp}]"), - ]: - ax.set_xlabel("h") - ax.set_ylabel(ylabel) - ax.set_title(title) - ax.legend() - ax.grid(True, alpha=0.3, which="both") - - fig_l2.tight_layout(); fig_h1.tight_layout() - fig_l2.savefig(os.path.join(results_dir, f"convergence_L2_{dim}_nu{nu}.png"), dpi=150) - fig_h1.savefig(os.path.join(results_dir, f"convergence_H1_{dim}_nu{nu}.png"), dpi=150) - plt.close(fig_l2); plt.close(fig_h1) - -# ============ MAIN ========================= -if __name__ == "__main__": - - print(f" MMS TYPE : {MMS_TYPE.upper()}") - - L = 1.0 - E = 1e6 - nx_values = [10, 20, 50, 80, 100, 120, 150, 200, 250 ] - - element_quad = _QuadElement() - element_tri = _TriElement() - - specs = [ - {"elem": element_quad, "label": "Q1 quad", "marker": "o", "color": "C0"}, - {"elem": element_tri, "label": "P1 tri", "marker": "s", "color": "C1"}, - ] - - # =========== Plane stress vec2d ======================= - DIM = "2d" - for nu in [0.0, 0.3, 0.49]: - print(f"\n PoissonRatio = {nu} ") - simulation_ponctuelle(element_quad, L, E, nu, nx=20, ny=20, dim=DIM) - simulation_ponctuelle(element_tri, L, E, nu, nx=20, ny=20, dim=DIM) - convergence_study(specs, L, E, nu, nx_values, dim=DIM) - - # ============== Plane strain vec3d ====================== - DIM = "3d" - - for nu in [0.0, 0.3, 0.49]: - print(f"\n PoissonRatio = {nu} ") - simulation_ponctuelle(element_quad, L, E, nu, nx=20, ny=20, dim=DIM) - simulation_ponctuelle(element_tri, L, E, nu, nx=20, ny=20, dim=DIM) - convergence_study(specs, L, E, nu, nx_values, dim=DIM) \ No newline at end of file diff --git a/examples/Freefem/MMS/2D/trigonometric.py b/examples/Freefem/MMS/2D/trigonometric.py new file mode 100644 index 00000000..6fea7bf8 --- /dev/null +++ b/examples/Freefem/MMS/2D/trigonometric.py @@ -0,0 +1,62 @@ +""" +Trigonometric 2D MMS on [0,L]^2 with linear-elasticity constitutive law: + + u_ex(x, y) = ( sin(pi x / L) cos(pi y / L), + cos(pi x / L) sin(pi y / L) ) + + sigma = lambda tr(eps) I + 2 mu eps, + with (lambda, mu) selected per dim (plane stress vs plane strain). +""" + +import numpy as np + +from manufactured_solution import MMSCase2D, lame + + +class Trigonometric(MMSCase2D): + name = "trigonometric" + plot_label = (r"$u_x = \sin(\pi x/L)\cos(\pi y/L),\ " + r"u_y = \cos(\pi x/L)\sin(\pi y/L)$") + + def u_ex(self, x, y, L): + k = np.pi / L + return (np.sin(k * x) * np.cos(k * y), + np.cos(k * x) * np.sin(k * y)) + + def grad_u_ex(self, x, y, L): + k = np.pi / L + dux_dx = k * np.cos(k * x) * np.cos(k * y) + dux_dy = -k * np.sin(k * x) * np.sin(k * y) + duy_dx = -k * np.sin(k * x) * np.sin(k * y) + duy_dy = k * np.cos(k * x) * np.cos(k * y) + return np.array([[dux_dx, dux_dy], + [duy_dx, duy_dy]]) + + def source(self, x, y, E, nu, L, dim): + lam, mu = lame(E, nu, dim) + k = np.pi / L + ux = np.sin(k * x) * np.cos(k * y) + uy = np.cos(k * x) * np.sin(k * y) + d2ux_dxx = -k**2 * ux + d2ux_dyy = -k**2 * ux + d2ux_dxy = -k**2 * np.cos(k * x) * np.sin(k * y) + d2uy_dxx = -k**2 * uy + d2uy_dyy = -k**2 * uy + d2uy_dxy = -k**2 * np.sin(k * x) * np.cos(k * y) + fx = -((lam + 2*mu) * d2ux_dxx + lam * d2uy_dxy + + mu * (d2ux_dyy + d2uy_dxy)) + fy = -(mu * (d2ux_dxy + d2uy_dxx) + lam * d2ux_dxy + + (lam + 2*mu) * d2uy_dyy) + return (fx, fy) + + +from beam import case_scene, element_quad, element_tri + +mms = Trigonometric() +createScene = case_scene(mms, element_quad) + + +if __name__ == "__main__": + from beam import run_reference_scene + + run_reference_scene(element_quad, mms) diff --git a/examples/Freefem/MMS/fem.py b/examples/Freefem/MMS/fem.py index 7a6c7482..146e2678 100644 --- a/examples/Freefem/MMS/fem.py +++ b/examples/Freefem/MMS/fem.py @@ -4,9 +4,6 @@ # --------------------------------------------------------------------------- # Quadrature rules -# Approximation: integral_{x_i}^{x_{i+1}} g(x) dx -# ~= (h/2) * sum_k( w_k * g(x_k) ) -# where x_k = (x_i + x_{i+1})/2 + (h/2)*xi_k # --------------------------------------------------------------------------- # Canonical Gauss-Legendre points and weights on the reference segment [-1, 1]. @@ -36,15 +33,106 @@ def rule(g, x1, x2): # Error-norm rules. H1 requires >=2 pts: the P1 gradient superconverges at # the element midpoint, so 1-pt fakes O(h^2) instead of O(h). -L2_QUADRATURE = line_quadrature(2) -H1_QUADRATURE = line_quadrature(2) +L2_QUADRATURE_1D = line_quadrature(2) +H1_QUADRATURE_1D = line_quadrature(2) + + +# --------------------------------------------------------------------------- +# 2D Q1 reference shape functions on [-1,1]^2. +# Node ordering: (-1,-1), (+1,-1), (+1,+1), (-1,+1) +# --------------------------------------------------------------------------- + +def _shape_q1(xi, eta): + 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)]) + 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 # --------------------------------------------------------------------------- # FEM assembly # --------------------------------------------------------------------------- -def assemble_nodal_forces(f_body, nodes, edges, quadrature): +def assemble_nodal_forces_1d(f_body, nodes, edges, quadrature): """ Assemble the consistent nodal force vector F_i = integral f_body(x) phi_i(x) dx. @@ -62,11 +150,72 @@ def assemble_nodal_forces(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 + + +def assemble_nodal_forces_2d(f_body, nodes, conn, element_rule): + """ + Assemble consistent nodal forces on a 2D mesh, agnostic to element type. + + F_a = sum_e integral_{Omega_e} N_a(x,y) * f_body(x,y) dx dy + + f_body : callable (x, y) -> (fx, fy) + nodes : (N, 2) or (N, 3) array — only the first two columns used + conn : iterable of node-index lists, one per element + element_rule : rule from quad_q1_rule(n_pts) / tri_p1_rule(n_pts) + """ + xy = np.asarray(nodes)[:, :2] + F = np.zeros((len(xy), 2)) + 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 + return F + + # --------------------------------------------------------------------------- # Error norms # --------------------------------------------------------------------------- -def l2_error(nodes, edges, u_h, u_ex, quadrature): +def l2_error_1d(nodes, edges, u_h, u_ex, quadrature): """L2 error norm: sqrt( integral (u_h - u_ex)^2 dx ) over the mesh.""" total = 0.0 for a, b in edges: @@ -78,7 +227,7 @@ def l2_error(nodes, edges, u_h, u_ex, quadrature): return np.sqrt(total) -def h1_semi_error(nodes, edges, u_h, du_ex, quadrature): +def h1_semi_error_1d(nodes, edges, u_h, du_ex, quadrature): """H1 semi-norm error: sqrt( integral (du_h - du_ex)^2 dx ) over the mesh.""" total = 0.0 for a, b in edges: @@ -88,3 +237,62 @@ def h1_semi_error(nodes, edges, u_h, du_ex, quadrature): du_h = u_h[a] * (-1.0 / h) + u_h[b] * (1.0 / h) total += quadrature(lambda x, du_h=du_h: (du_h - du_ex(x)) ** 2, x1, x2) 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) diff --git a/examples/Freefem/MMS/plot_utils.py b/examples/Freefem/MMS/plot_utils.py new file mode 100644 index 00000000..5764c32d --- /dev/null +++ b/examples/Freefem/MMS/plot_utils.py @@ -0,0 +1,16 @@ +"""Plot helpers shared across MMS 1D/2D drivers.""" + +import numpy as np + + +def annotate_convergence_rates(ax, hs, errors): + """Annotate the per-segment log-log convergence rate above each segment.""" + hs = np.asarray(hs) + errors = np.asarray(errors) + for k in range(1, len(hs)): + rate = np.log(errors[k] / errors[k-1]) / np.log(hs[k] / hs[k-1]) + x_mid = np.sqrt(hs[k] * hs[k-1]) + y_mid = np.sqrt(errors[k] * errors[k-1]) + ax.annotate(f"{rate:.2f}", xy=(x_mid, y_mid), + xytext=(0, 8), textcoords="offset points", + ha="center", fontsize=9) From fc7353710700e5f78a1134efcac3aee39aa7ab99 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 26 May 2026 15:19:39 +0200 Subject: [PATCH 2/4] Factorize BCs into function scripts --- examples/Freefem/MMS/1D/bar.py | 6 +--- examples/Freefem/MMS/1D/cubic.py | 7 ++++ examples/Freefem/MMS/1D/exponential.py | 7 ++++ .../Freefem/MMS/1D/manufactured_solution.py | 9 +++++ examples/Freefem/MMS/1D/quadratic.py | 7 ++++ examples/Freefem/MMS/1D/run_convergence.py | 17 ++------- examples/Freefem/MMS/1D/sinusoidal.py | 7 ++++ examples/Freefem/MMS/2D/beam.py | 36 +------------------ examples/Freefem/MMS/2D/cubic.py | 31 ++++++++++++++++ .../Freefem/MMS/2D/manufactured_solution.py | 10 ++++++ examples/Freefem/MMS/2D/run_convergence.py | 24 ++----------- examples/Freefem/MMS/2D/trigonometric.py | 31 ++++++++++++++++ 12 files changed, 116 insertions(+), 76 deletions(-) diff --git a/examples/Freefem/MMS/1D/bar.py b/examples/Freefem/MMS/1D/bar.py index 8a6c2f41..b1fe7b50 100644 --- a/examples/Freefem/MMS/1D/bar.py +++ b/examples/Freefem/MMS/1D/bar.py @@ -130,11 +130,7 @@ 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): diff --git a/examples/Freefem/MMS/1D/cubic.py b/examples/Freefem/MMS/1D/cubic.py index eb025550..0ad96e67 100644 --- a/examples/Freefem/MMS/1D/cubic.py +++ b/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) diff --git a/examples/Freefem/MMS/1D/exponential.py b/examples/Freefem/MMS/1D/exponential.py index d5c50a8c..8f0d74e1 100644 --- a/examples/Freefem/MMS/1D/exponential.py +++ b/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) diff --git a/examples/Freefem/MMS/1D/manufactured_solution.py b/examples/Freefem/MMS/1D/manufactured_solution.py index 058a83eb..1f8c9789 100644 --- a/examples/Freefem/MMS/1D/manufactured_solution.py +++ b/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)`. + """ diff --git a/examples/Freefem/MMS/1D/quadratic.py b/examples/Freefem/MMS/1D/quadratic.py index a0c2fb6a..85a8c9fa 100644 --- a/examples/Freefem/MMS/1D/quadratic.py +++ b/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) diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py index b68615af..0e3a51f0 100644 --- a/examples/Freefem/MMS/1D/run_convergence.py +++ b/examples/Freefem/MMS/1D/run_convergence.py @@ -87,35 +87,22 @@ def write_convergence_table(stem, rows): f.write(line + "\n") -def plot_convergence(stem, hs, series, title, ylabel="Error", ref_lines=None): +def plot_convergence(stem, hs, series, title, ylabel="Error"): """ Save log-log convergence plot to results/.png. - series : list of {"label", "errors", "style"?} dicts - ref_lines : optional list of {"slope", "anchor", "label"?} dicts; - "anchor" names the series whose first point anchors the line. + 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)) - anchors = {} 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) - anchors[s["label"]] = e_arr - for rl in (ref_lines or []): - anchor = anchors.get(rl["anchor"]) - if anchor is None: - continue - slope = rl["slope"] - label = rl.get("label", f"O(h^{slope})") - h_ref = np.array([h_arr[0], h_arr[-1]]) - ax.loglog(h_ref, anchor[0] * (h_ref / h_arr[0])**slope, ":", - color="gray", lw=1.2, label=label) ax.set_xlabel("h") ax.set_ylabel(ylabel) ax.set_title(title) diff --git a/examples/Freefem/MMS/1D/sinusoidal.py b/examples/Freefem/MMS/1D/sinusoidal.py index 317e7f71..63805a5b 100644 --- a/examples/Freefem/MMS/1D/sinusoidal.py +++ b/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) diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py index 0e2046db..4fe58f82 100644 --- a/examples/Freefem/MMS/2D/beam.py +++ b/examples/Freefem/MMS/2D/beam.py @@ -52,40 +52,6 @@ def _dim_template(dim): return "Vec3d" if dim == "3d" else "Vec2d" -def _add_dirichlet(Beam, nodes_2d, L, dim, with_visual): - """Partial Dirichlet: ux=0 on x-faces, uy=0 on y-faces (and z=0 in 3d).""" - tmpl = _dim_template(dim) - xy = nodes_2d[:, :2] - eps = 1e-10 - - idx_ux0 = [k for k, (xk, _) in enumerate(xy) if xk < eps or xk > L - eps] - idx_uy0 = [k for k, (_, yk) in enumerate(xy) if yk < eps or yk > L - eps] - - if dim == "2d": - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_ux", template=tmpl, - indices=" ".join(map(str, idx_ux0)), - fixedDirections="1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_uy", template=tmpl, - indices=" ".join(map(str, idx_uy0)), - fixedDirections="0 1") - else: - all_idx = " ".join(map(str, range(len(nodes_2d)))) - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_ux", template=tmpl, - indices=" ".join(map(str, idx_ux0)), - fixedDirections="1 0 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_uy", template=tmpl, - indices=" ".join(map(str, idx_uy0)), - fixedDirections="0 1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_z", template=tmpl, - indices=all_idx, - fixedDirections="0 0 1") - - 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)] @@ -300,7 +266,7 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, Beam.addObject("LinearSmallStrainFEMForceField", name="FEM", template=tmpl, youngModulus=E, poissonRatio=nu, topology="@topology") - _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) + mms.apply_bcs(Beam, nodes_2d, L, dim) # Placeholder force field filled in by the controller after init. n_nodes = len(nodes_2d) diff --git a/examples/Freefem/MMS/2D/cubic.py b/examples/Freefem/MMS/2D/cubic.py index 43c17afc..7c94168b 100644 --- a/examples/Freefem/MMS/2D/cubic.py +++ b/examples/Freefem/MMS/2D/cubic.py @@ -29,6 +29,37 @@ def grad_u_ex(self, x, y, L): return np.array([[dux_dx, dux_dy], [duy_dx, duy_dy]]) + def apply_bcs(self, Beam, nodes_2d, L, dim): + tmpl = "Vec3d" if dim == "3d" else "Vec2d" + xy = nodes_2d[:, :2] + eps = 1e-10 + idx_ux0 = [k for k, (xk, _) in enumerate(xy) if xk < eps or xk > L - eps] + idx_uy0 = [k for k, (_, yk) in enumerate(xy) if yk < eps or yk > L - eps] + + if dim == "2d": + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1") + else: + all_idx = " ".join(map(str, range(len(nodes_2d)))) + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_z", template=tmpl, + indices=all_idx, + fixedDirections="0 0 1") + def source(self, x, y, E, nu, L, dim): lam, mu = lame(E, nu, dim) d2ux_dxx = 2*L - 6*x - 2*y diff --git a/examples/Freefem/MMS/2D/manufactured_solution.py b/examples/Freefem/MMS/2D/manufactured_solution.py index 88ce9408..573e9bce 100644 --- a/examples/Freefem/MMS/2D/manufactured_solution.py +++ b/examples/Freefem/MMS/2D/manufactured_solution.py @@ -39,6 +39,16 @@ def source(self, x, y, E, nu, L, dim): "2d" → plane stress, "3d" → plane strain. """ + @abstractmethod + def apply_bcs(self, Beam, nodes_2d, L, dim): + """Install Dirichlet BCs (and any case-specific constraints) on `Beam`. + + The MMS author chooses the BC pattern matching the manufactured + solution. Neumann tractions on the four edges are already assembled + into the consistent nodal force by the framework (via `self.traction`); + this method only needs to add SOFA constraint objects. + """ + def traction(self, x, y, nx, ny, E, nu, L, dim): """Traction on a face with outward unit normal (nx, ny): returns (Tx, Ty). diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py index 724b73d3..ad8ad0a5 100644 --- a/examples/Freefem/MMS/2D/run_convergence.py +++ b/examples/Freefem/MMS/2D/run_convergence.py @@ -70,14 +70,9 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"): "errors": h1s, "style": spec["h1_style"]}) hs_ref = hs - first_label = elem_specs[0]["label"] - ref_lines = [ - {"slope": 2, "anchor": f"{first_label} L²", "label": "O(h²)"}, - {"slope": 1, "anchor": f"{first_label} H¹", "label": "O(h¹)"}, - ] title = f"Convergence — {mms.name} [{dim} / {hyp}] nu={nu}" plot_convergence(f"convergence_{mms.name}_{dim}_nu{nu}", - hs_ref, plot_series, title=title, ref_lines=ref_lines) + hs_ref, plot_series, title=title) def write_convergence_table(stem, rows): @@ -102,35 +97,22 @@ def write_convergence_table(stem, rows): f.write(line + "\n") -def plot_convergence(stem, hs, series, title, ylabel="Error", ref_lines=None): +def plot_convergence(stem, hs, series, title, ylabel="Error"): """ Save log-log convergence plot to results/.png. - series : list of {"label", "errors", "style"?} dicts - ref_lines : optional list of {"slope", "anchor", "label"?} dicts; - "anchor" names the series whose first point anchors the line. + 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)) - anchors = {} 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) - anchors[s["label"]] = e_arr - for rl in (ref_lines or []): - anchor = anchors.get(rl["anchor"]) - if anchor is None: - continue - slope = rl["slope"] - label = rl.get("label", f"O(h^{slope})") - h_ref = np.array([h_arr[0], h_arr[-1]]) - ax.loglog(h_ref, anchor[0] * (h_ref / h_arr[0])**slope, ":", - color="gray", lw=1.2, label=label) ax.set_xlabel("h") ax.set_ylabel(ylabel) ax.set_title(title) diff --git a/examples/Freefem/MMS/2D/trigonometric.py b/examples/Freefem/MMS/2D/trigonometric.py index 6fea7bf8..e9527a33 100644 --- a/examples/Freefem/MMS/2D/trigonometric.py +++ b/examples/Freefem/MMS/2D/trigonometric.py @@ -32,6 +32,37 @@ def grad_u_ex(self, x, y, L): return np.array([[dux_dx, dux_dy], [duy_dx, duy_dy]]) + def apply_bcs(self, Beam, nodes_2d, L, dim): + tmpl = "Vec3d" if dim == "3d" else "Vec2d" + xy = nodes_2d[:, :2] + eps = 1e-10 + idx_ux0 = [k for k, (xk, _) in enumerate(xy) if xk < eps or xk > L - eps] + idx_uy0 = [k for k, (_, yk) in enumerate(xy) if yk < eps or yk > L - eps] + + if dim == "2d": + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1") + else: + all_idx = " ".join(map(str, range(len(nodes_2d)))) + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_ux", template=tmpl, + indices=" ".join(map(str, idx_ux0)), + fixedDirections="1 0 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_uy", template=tmpl, + indices=" ".join(map(str, idx_uy0)), + fixedDirections="0 1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_z", template=tmpl, + indices=all_idx, + fixedDirections="0 0 1") + def source(self, x, y, E, nu, L, dim): lam, mu = lame(E, nu, dim) k = np.pi / L From 07dc5032a56bdfc87eee5e1fafaa29e5962fa487 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 26 May 2026 15:22:34 +0200 Subject: [PATCH 3/4] Pull incompressible case under the hood --- examples/Freefem/2D/neumann_incompressible.py | 644 ------------------ examples/Freefem/MMS/2D/incompressible.py | 102 +++ examples/Freefem/MMS/2D/params.json | 7 +- examples/Freefem/MMS/2D/run_convergence.py | 7 +- 4 files changed, 110 insertions(+), 650 deletions(-) delete mode 100644 examples/Freefem/2D/neumann_incompressible.py create mode 100644 examples/Freefem/MMS/2D/incompressible.py diff --git a/examples/Freefem/2D/neumann_incompressible.py b/examples/Freefem/2D/neumann_incompressible.py deleted file mode 100644 index e3c43cfc..00000000 --- a/examples/Freefem/2D/neumann_incompressible.py +++ /dev/null @@ -1,644 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import os - -import Sofa -import Sofa.Core -import Sofa.Simulation -import SofaRuntime - -""" -Incompressible ===> div u(x,y) = 0, -MMS realistic : u_x = sin(pi*x)cos(pi*y); u_y = -cos(pi*x)sin(pi*y) -""" - -RESULTS_DIR = "results_incompressible_q1p1" -os.makedirs(RESULTS_DIR, exist_ok=True) - -# ===== Gauss points & weights ===== -GAUSS_PTS = np.array([-1.0 / np.sqrt(3), 1.0 / np.sqrt(3)]) -GAUSS_WTS = np.array([1.0, 1.0]) - -# ===== Triangle Gauss rule ===== -_TRI_BARY = np.array([ - [2/3, 1/6, 1/6], - [1/6, 2/3, 1/6], - [1/6, 1/6, 2/3], -]) - - -# ======= Lamé coefficients ========================= -def lame(E, nu, dim="2d"): - if dim == "3d" and np.isclose(nu, 0.5): - raise ValueError( - "nu=0.5 in plane strain ===> division by zero (strictly incompressible). " - "Use nu < 0.5 " - ) - if dim == "2d": - lam = E * nu / (1 - nu**2) - else: - lam = E * nu / ((1 + nu) * (1 - 2 * nu)) - mu = E / (2 * (1 + nu)) - return lam, mu - - -# ==================== MMS & derivatives ========================== -def ux_mms_incomp(x, y, L): return np.sin(np.pi*x/L) * np.cos(np.pi*y/L) -def uy_mms_incomp(x, y, L): return -np.cos(np.pi*x/L) * np.sin(np.pi*y/L) - -def dux_dx_incomp(x, y, L): return (np.pi/L)*np.cos(np.pi*x/L)*np.cos(np.pi*y/L) -def dux_dy_incomp(x, y, L): return -(np.pi/L)*np.sin(np.pi*x/L)*np.sin(np.pi*y/L) -def duy_dx_incomp(x, y, L): return (np.pi/L)*np.sin(np.pi*x/L)*np.sin(np.pi*y/L) -def duy_dy_incomp(x, y, L): return -(np.pi/L)*np.cos(np.pi*x/L)*np.cos(np.pi*y/L) - -def fx_body_incomp(x, y, E, nu, L, dim="2d"): - return (np.pi**2 * E / (L**2 * (nu + 1))) * np.sin(np.pi*x/L) * np.cos(np.pi*y/L) - -def fy_body_incomp(x, y, E, nu, L, dim="2d"): - return -(np.pi**2 * E / (L**2 * (nu + 1))) * np.sin(np.pi*y/L) * np.cos(np.pi*x/L) - -def sigma_xx_incomp(x, y, E, nu, L, dim="2d"): - _, mu = lame(E, nu, dim) - return 2*mu * dux_dx_incomp(x, y, L) - -def sigma_yy_incomp(x, y, E, nu, L, dim="2d"): - _, mu = lame(E, nu, dim) - return 2*mu * duy_dy_incomp(x, y, L) - -def sigma_xy_incomp(x, y, E, nu, L, dim="2d"): - _, mu = lame(E, nu, dim) - return mu * (dux_dy_incomp(x, y, L) + duy_dx_incomp(x, y, L)) - -def traction_incomp(x, y, nx_c, ny_c, E, nu, L, dim="2d"): - sxx = sigma_xx_incomp(x, y, E, nu, L, dim) - syy = sigma_yy_incomp(x, y, E, nu, L, dim) - sxy = sigma_xy_incomp(x, y, E, nu, L, dim) - return sxx*nx_c + sxy*ny_c, sxy*nx_c + syy*ny_c - - -# =================== Mesh ============================================ -def get_nodes_2d(L, nx, ny, dim="2d"): - dx, dy = L / (nx - 1), L / (ny - 1) - if dim == "3d": - pts = [[i*dx, j*dy, 0.0] for j in range(ny) for i in range(nx)] - else: - pts = [[i*dx, j*dy] for j in range(ny) for i in range(nx)] - return np.array(pts, dtype=float) - -def _sofa_template(dim): - return "Vec3d" if dim == "3d" else "Vec2d" - -def _add_dirichlet(Beam, nodes_2d, L, dim, with_visual): - tmpl = _sofa_template(dim) - xy = nodes_2d[:, :2] - eps = 1e-10 - - idx_corner = [k for k, (xk, yk) in enumerate(xy) - if xk < eps and yk < eps] - idx_rot = [k for k, (xk, yk) in enumerate(xy) - if xk > L - eps and yk < eps] - - if dim == "2d": - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_corner", template=tmpl, - indices=" ".join(map(str, idx_corner)), - fixedDirections="1 1") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_rot", template=tmpl, - indices=" ".join(map(str, idx_rot)), - fixedDirections="0 1") - else: - all_idx = " ".join(str(k) for k in range(len(nodes_2d))) - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_corner", template=tmpl, - indices=" ".join(map(str, idx_corner)), - fixedDirections="1 1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_rot", template=tmpl, - indices=" ".join(map(str, idx_rot)), - fixedDirections="0 1 0") - Beam.addObject("PartialFixedProjectiveConstraint", - name="fix_z", template=tmpl, - indices=all_idx, - fixedDirections="0 0 1") - -def _pack_forces(F_xy, dim, n_nodes): - idx = " ".join(str(k) for k in range(n_nodes)) - if dim == "3d": - frc = " ".join(f"{F_xy[k,0]} {F_xy[k,1]} 0.0" for k in range(n_nodes)) - else: - frc = " ".join(f"{F_xy[k,0]} {F_xy[k,1]}" for k in range(n_nodes)) - return idx, frc - - -# ======================== Q1 Quad element ============================================= - -class _QuadElement: - LABEL = "Q1 quad" - - @staticmethod - def get_connectivity(nx, ny): - quads = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - quads.append([k00, k10, k11, k01]) - return np.array(quads) - - @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) - - @staticmethod - def _shape(xi, eta): - 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)]) - return N, dN_dxi, dN_deta - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim="2d"): - quads = _QuadElement.get_connectivity(nx, ny) - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - - - for quad in quads: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(xi, eta) - J = np.array([[dN_dxi@xe, dN_dxi@ye], - [dN_deta@xe, dN_deta@ye]]) - detJ = np.linalg.det(J) - xg, yg = N@xe, N@ye - fx = fx_body_incomp(xg, yg, E, nu, L, dim) - fy = fy_body_incomp(xg, yg, E, nu, L, dim) - w = wi * wj * detJ - for a, node in enumerate(quad): - F[node, 0] += N[a] * fx * w - F[node, 1] += N[a] * fy * w - - - def _edge_gauss(k0, k1, nx_c, ny_c): - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5*(xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction_incomp(xg, yg, nx_c, ny_c, E, nu, L, dim) - N0, N1 = 0.5*(1-xi), 0.5*(1+xi) - w = wi * Le / 2.0 - F[k0, 0] += N0*Tx*w; F[k0, 1] += N0*Ty*w - F[k1, 0] += N1*Tx*w; F[k1, 1] += N1*Ty*w - - for j in range(ny - 1): # x=0 - _edge_gauss(j*nx, (j+1)*nx, -1., 0.) - for j in range(ny - 1): # x=L - _edge_gauss(j*nx+(nx-1), (j+1)*nx+(nx-1), +1., 0.) - for i in range(nx - 1): # y=0 - _edge_gauss(i, i+1, 0., -1.) - for i in range(nx - 1): # y=L - _edge_gauss((ny-1)*nx+i, (ny-1)*nx+i+1, 0., +1.) - - return F - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nu=0.3, - nx=10, ny=10, with_visual=True, dim="2d"): - tmpl = _sofa_template(dim) - rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - quads_np = _QuadElement.get_connectivity(nx, ny) - - # ========== SOFA solvers ====================== - - - Beam = rootNode.addChild("Beam") - Beam.addObject('NewtonRaphsonSolver', - name="newtonSolver", - maxNbIterationsNewton=10, - absoluteResidualStoppingThreshold=1e-10, - printLog=False) - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), - showObject=with_visual, showObjectScale=0.005*L) - Beam.addObject("QuadSetTopologyContainer", - name="topology", quads=quads_np.tolist()) - Beam.addObject("QuadSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", - name="FEM", template=tmpl, - youngModulus=E, poissonRatio=nu, topology="@topology") - - _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) - - F_xy = _QuadElement.compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim) - idx, frc = _pack_forces(F_xy, dim, len(nodes_2d)) - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - return dofs, nodes_2d, quads_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, conn): - xy, err2 = nodes_2d[:, :2], 0.0 - for quad in conn: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(xi, eta) - J = np.array([[dN_dxi@xe, dN_dxi@ye], - [dN_deta@xe, dN_deta@ye]]) - detJ = np.linalg.det(J) - xg, yg = N@xe, N@ye - err2 += ( - (N@ux[quad] - ux_mms_incomp(xg, yg, L))**2 - + (N@uy[quad] - uy_mms_incomp(xg, yg, L))**2 - ) * wi * wj * detJ - return np.sqrt(err2) - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, conn, **_): - xy, err2 = nodes_2d[:, :2], 0.0 - for quad in conn: - xe, ye = xy[quad, 0], xy[quad, 1] - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - for eta, wj in zip(GAUSS_PTS, GAUSS_WTS): - N, dN_dxi, dN_deta = _QuadElement._shape(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 - err2 += ( - (dN_dx@ux[quad] - dux_dx_incomp(xg, yg, L))**2 + - (dN_dy@ux[quad] - dux_dy_incomp(xg, yg, L))**2 + - (dN_dx@uy[quad] - duy_dx_incomp(xg, yg, L))**2 + - (dN_dy@uy[quad] - duy_dy_incomp(xg, yg, L))**2 - ) * wi * wj * detJ - return np.sqrt(err2) - - - -# ==================== P1 Triangle element ========================== - -class _TriElement: - LABEL = "P1 tri" - - @staticmethod - def get_connectivity(nx, ny): - tris = [] - for j in range(ny - 1): - for i in range(nx - 1): - k00 = j * nx + i - k10 = j * nx + (i + 1) - k01 = (j + 1) * nx + i - k11 = (j + 1) * nx + (i + 1) - tris += [[k00, k10, k11], [k00, k11, k01]] - return np.array(tris) - - @staticmethod - def to_triangles(conn): - return conn - - @staticmethod - def _grad_p1(xy, u, tri): - i0, i1, i2 = tri - x0,y0 = xy[i0]; x1,y1 = xy[i1]; x2,y2 = xy[i2] - A2 = (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0) - dudx = (u[i0]*(y1-y2) + u[i1]*(y2-y0) + u[i2]*(y0-y1)) / A2 - dudy = (u[i0]*(x2-x1) + u[i1]*(x0-x2) + u[i2]*(x1-x0)) / A2 - return dudx, dudy, abs(A2) / 2.0 - - @staticmethod - def compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim="2d"): - tris = _TriElement.get_connectivity(nx, ny) - xy = nodes_2d[:, :2] - F = np.zeros((len(xy), 2)) - - for tri in tris: - i0, i1, i2 = tri - x0,y0 = xy[i0]; x1,y1 = xy[i1]; x2,y2 = xy[i2] - xc = (x0 + x1 + x2) / 3.0 - yc = (y0 + y1 + y2) / 3.0 - area = abs((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)) / 2.0 - fx = fx_body_incomp(xc, yc, E, nu, L, dim) - fy = fy_body_incomp(xc, yc, E, nu, L, dim) - for node in [i0, i1, i2]: - F[node, 0] += fx * area / 3.0 - F[node, 1] += fy * area / 3.0 - - - def _edge_gauss(k0, k1, nx_c, ny_c): - x0, y0 = xy[k0]; x1, y1 = xy[k1] - Le = np.hypot(x1-x0, y1-y0) - for xi, wi in zip(GAUSS_PTS, GAUSS_WTS): - t = 0.5*(xi + 1) - xg, yg = (1-t)*x0 + t*x1, (1-t)*y0 + t*y1 - Tx, Ty = traction_incomp(xg, yg, nx_c, ny_c, E, nu, L, dim) - N0, N1 = 0.5*(1-xi), 0.5*(1+xi) - w = wi * Le / 2.0 - F[k0, 0] += N0*Tx*w; F[k0, 1] += N0*Ty*w - F[k1, 0] += N1*Tx*w; F[k1, 1] += N1*Ty*w - - for j in range(ny - 1): - _edge_gauss(j*nx, (j+1)*nx, -1., 0.) - for j in range(ny - 1): - _edge_gauss(j*nx+(nx-1), (j+1)*nx+(nx-1), +1., 0.) - for i in range(nx - 1): - _edge_gauss(i, i+1, 0., -1.) - for i in range(nx - 1): - _edge_gauss((ny-1)*nx+i, (ny-1)*nx+i+1, 0., +1.) - - return F - - @staticmethod - def create_sofa_scene(rootNode, L=1.0, E=1e6, nu=0.3, - nx=10, ny=10, with_visual=True, dim="2d"): - tmpl = _sofa_template(dim) - rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") - rootNode.addObject("RequiredPlugin", pluginName=[ - "Elasticity", "Sofa.Component.Constraint.Projective", - "Sofa.Component.Engine.Select", "Sofa.Component.LinearSolver.Direct", - "Sofa.Component.MechanicalLoad", "Sofa.Component.ODESolver.Backward", - "Sofa.Component.StateContainer", - "Sofa.Component.Topology.Container.Dynamic", - ]) - rootNode.addObject("DefaultAnimationLoop") - if with_visual: - rootNode.addObject("VisualStyle", - displayFlags="showBehaviorModels showForceFields") - - nodes_2d = get_nodes_2d(L, nx, ny, dim=dim) - tris_np = _TriElement.get_connectivity(nx, ny) - - Beam = rootNode.addChild("Beam") - Beam.addObject("StaticSolver", name="staticSolver", printLog=False) - Beam.addObject("SparseLDLSolver", name="linearSolver", - template="CompressedRowSparseMatrixd") - dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl, - position=nodes_2d.tolist(), - showObject=with_visual, showObjectScale=0.005*L) - Beam.addObject("TriangleSetTopologyContainer", - name="topology", triangles=tris_np.tolist()) - Beam.addObject("TriangleSetTopologyModifier") - Beam.addObject("LinearSmallStrainFEMForceField", - name="FEM", template=tmpl, - youngModulus=E, poissonRatio=nu, topology="@topology") - - _add_dirichlet(Beam, nodes_2d, L, dim, with_visual) - - F_xy = _TriElement.compute_nodal_forces(nodes_2d, L, E, nu, nx, ny, dim) - idx, frc = _pack_forces(F_xy, dim, len(nodes_2d)) - Beam.addObject("ConstantForceField", name="MMS_forces", template=tmpl, - indices=idx, forces=frc) - return dofs, nodes_2d, tris_np - - @staticmethod - def compute_l2(nodes_2d, ux, uy, L, nx, ny, conn): - - xy = nodes_2d[:, :2] - err2 = 0.0 - for tri in conn: - i0, i1, i2 = tri - xy_e = xy[[i0, i1, i2]] - ux_e = ux[[i0, i1, i2]] - uy_e = uy[[i0, i1, i2]] - _, _, area = _TriElement._grad_p1(xy, ux, tri) - for lam in _TRI_BARY: - xg = lam @ xy_e[:, 0] - yg = lam @ xy_e[:, 1] - uxg = lam @ ux_e - uyg = lam @ uy_e - err2 += ( - (uxg - ux_mms_incomp(xg, yg, L))**2 - + (uyg - uy_mms_incomp(xg, yg, L))**2 - ) * area / 3.0 - return np.sqrt(err2) - - @staticmethod - def compute_h1(nodes_2d, ux, uy, L, conn, **_): - xy = nodes_2d[:, :2] - err2 = 0.0 - for tri in conn: - i0, i1, i2 = tri - xc = (xy[i0,0] + xy[i1,0] + xy[i2,0]) / 3.0 - yc = (xy[i0,1] + xy[i1,1] + xy[i2,1]) / 3.0 - dux_dx_h, dux_dy_h, area = _TriElement._grad_p1(xy, ux, tri) - duy_dx_h, duy_dy_h, _ = _TriElement._grad_p1(xy, uy, tri) - err2 += ( - (dux_dx_h - dux_dx_incomp(xc, yc, L))**2 + - (dux_dy_h - dux_dy_incomp(xc, yc, L))**2 + - (duy_dx_h - duy_dx_incomp(xc, yc, L))**2 + - (duy_dy_h - duy_dy_incomp(xc, yc, L))**2 - ) * area - return np.sqrt(err2) - - -# ============================== SIMU======================================= - -def run_simulation(elem, L, E, nu, nx, ny, dim="2d"): - root = Sofa.Core.Node("root") - dofs, nodes_2d, conn = elem.create_sofa_scene( - root, L=L, E=E, nu=nu, nx=nx, ny=ny, with_visual=False, dim=dim - ) - Sofa.Simulation.init(root) - pos0 = dofs.position.array().copy() - Sofa.Simulation.animate(root, root.dt.value) - pos1 = dofs.position.array().copy() - Sofa.Simulation.unload(root) - ux = pos1[:, 0] - pos0[:, 0] - uy = pos1[:, 1] - pos0[:, 1] - return nodes_2d, conn, ux, uy - - -def simulation_ponctuelle(elem, L, E, nu, nx, ny, dim="2d", results_dir=RESULTS_DIR): - os.makedirs(results_dir, exist_ok=True) - nodes_2d, conn, ux, uy = run_simulation(elem, L, E, nu, nx, ny, dim=dim) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - - xy = nodes_2d[:, :2] - ux_ref = ux_mms_incomp(xy[:, 0], xy[:, 1], L) - uy_ref = uy_mms_incomp(xy[:, 0], xy[:, 1], L) - - hyp = "plane strain" if dim == "3d" else "plane stress" - print(f" L2 = {l2:.4e}") - print(f" H1 semi-norm = {h1:.4e}") - print(f" max|ux-ux_mms| = {np.max(np.abs(ux - ux_ref)):.4e}") - print(f" max|uy-uy_mms| = {np.max(np.abs(uy - uy_ref)):.4e}") - - mid_j = (ny - 1) // 2 - sl = slice(mid_j * nx, mid_j * nx + nx) - yc = xy[mid_j * nx, 1] - xf = np.linspace(0, L, 300) - - fig, axes = plt.subplots(1, 2, figsize=(12, 4)) - for ax, u_sofa, u_fn, lbl, fmt in zip( - axes, - [ux[sl], uy[sl]], - [lambda x: ux_mms_incomp(x, yc, L), lambda x: uy_mms_incomp(x, yc, L)], - [r"$u_x$", r"$u_y$"], ["o-", "s-"], - ): - ax.plot(xy[sl, 0], u_sofa, fmt, color="tab:orange", - label=f"SOFA {elem.LABEL} [{dim}]", ms=5) - ax.plot(xf, u_fn(xf), "--", color="tab:green", label="MMS exact") - ax.set_xlabel("x"); ax.set_ylabel(lbl) - ax.legend(); ax.grid(True, alpha=0.3) - plt.suptitle( - f"MMS INCOMPRESSIBLE — {elem.LABEL} [{dim}/{hyp}] nu={nu} nx={nx}" - f" |L2={l2:.2e} H1={h1:.2e}" - ) - plt.tight_layout() - tag = elem.LABEL.replace(" ", "_") - plt.savefig(os.path.join(results_dir, - f"solution_{tag}_{dim}_nu{nu}_nx{nx}.png"), dpi=150) - plt.close() - - tris_plot = elem.to_triangles(conn) - x, y = xy[:, 0], xy[:, 1] - fig, axes = plt.subplots(2, 3, figsize=(15, 8)) - for ax, data, title, cmap in [ - (axes[0,0], ux, r"$u_x$ SOFA", "gray"), - (axes[0,1], ux_ref, r"$u_x$ MMS", "gray"), - (axes[0,2], np.abs(ux-ux_ref), r"$|u_x-u_x^{MMS}|$", "gray_r"), - (axes[1,0], uy, r"$u_y$ SOFA", "gray"), - (axes[1,1], uy_ref, r"$u_y$ MMS", "gray"), - (axes[1,2], np.abs(uy-uy_ref), r"$|u_y-u_y^{MMS}|$", "gray_r"), - ]: - tc = ax.tricontourf(x, y, tris_plot.tolist(), data, levels=20, cmap=cmap) - ax.triplot(x, y, tris_plot.tolist(), "k-", lw=0.3, alpha=0.4) - plt.colorbar(tc, ax=ax, shrink=0.8) - ax.set_title(title); ax.set_aspect("equal") - ax.set_xlabel("x"); ax.set_ylabel("y") - plt.suptitle( - f"Champs 2D — {elem.LABEL} [{dim}/{hyp}] INCOMPRESSIBLE nu={nu} nx={nx}" - ) - plt.tight_layout() - plt.savefig(os.path.join(results_dir, - f"champs2D_{tag}_{dim}_nu{nu}_nx{nx}.png"), dpi=150) - plt.close() - - -def convergence_study(elem_specs, L, E, nu, nx_values, dim="2d", results_dir=RESULTS_DIR): - os.makedirs(results_dir, exist_ok=True) - hyp = "plane strain" if dim == "3d" else "plane stress" - fig_l2, ax_l2 = plt.subplots(figsize=(7, 5)) - fig_h1, ax_h1 = plt.subplots(figsize=(7, 5)) - hs_last = l2_last = h1_last = None - - for spec in elem_specs: - elem = spec["elem"] - label = spec["label"] - marker = spec.get("marker", "o") - color = spec.get("color", None) - - hs, errs_l2, errs_h1 = [], [], [] - hdr = (f"{'nx':>5} | {'h':>8} | {'L2':>14} | {'ord_L2':>7}" - f" | {'H1':>14} | {'ord_H1':>7}") - tag = label.replace(" ", "_") - txt_path = os.path.join(results_dir, - f"convergence_{tag}_{dim}_nu{nu}.txt") - - print(f"\n── Convergence {label} [{dim}/{hyp}] INCOMP nu={nu} ──\n{hdr}") - with open(txt_path, "w", encoding="utf-8") as f: - f.write(f"Convergence MMS INCOMPRESSIBLE — {label} [{dim}/{hyp}] nu={nu}\n{hdr}\n") - for k, nx in enumerate(nx_values): - ny = nx - h = L / (nx - 1) - nodes_2d, conn, ux, uy = run_simulation(elem, L, E, nu, nx, ny, dim=dim) - l2 = elem.compute_l2(nodes_2d, ux, uy, L, nx, ny, conn) - h1 = elem.compute_h1(nodes_2d, ux, uy, L, conn) - hs.append(h); errs_l2.append(l2); errs_h1.append(h1) - - ord_l2 = (f"{np.log(l2/errs_l2[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else " — ") - ord_h1 = (f"{np.log(h1/errs_h1[k-1])/np.log(h/hs[k-1]):.2f}" - if k > 0 else " — ") - line = (f"{nx:5d} | {h:8.4f} | {l2:14.6e} | {ord_l2:>7}" - f" | {h1:14.6e} | {ord_h1:>7}") - print(line); f.write(line + "\n") - - hs_a = np.array(hs); l2_a = np.array(errs_l2); h1_a = np.array(errs_h1) - kw = dict(lw=2, ms=7, color=color) - ax_l2.loglog(hs_a, l2_a, f"{marker}-", - label=f"{label} [{dim}] INCOMPRESSIBLE nu={nu}", **kw) - ax_h1.loglog(hs_a, h1_a, f"{marker}--", - label=f"{label} [{dim}] INCOMPRESSIBLE nu={nu}", **kw) - hs_last, l2_last, h1_last = hs_a, l2_a, h1_a - - if hs_last is not None: - h_ref = np.array([hs_last[0], hs_last[-1]]) - ax_l2.loglog(h_ref, l2_last[0]*(h_ref/hs_last[0])**2, - ":", color="gray", lw=1.2, label="O(h²)") - ax_h1.loglog(h_ref, h1_last[0]*(h_ref/hs_last[0])**1, - ":", color="gray", lw=1.2, label="O(h¹)") - - for ax, ylabel, title in [ - (ax_l2, "L2 error", - f"Convergence L2 — MMS INCOMPRESSIBLE [{dim}/{hyp}]"), - (ax_h1, "H1 semi-norm", - f"Convergence H1 — MMS INCOMPRESSIBLE [{dim}/{hyp}]"), - ]: - ax.set_xlabel("h"); ax.set_ylabel(ylabel); ax.set_title(title) - ax.legend(); ax.grid(True, alpha=0.3, which="both") - - fig_l2.tight_layout(); fig_h1.tight_layout() - fig_l2.savefig(os.path.join(results_dir, - f"convergence_L2_{dim}_nu{nu}.png"), dpi=150) - fig_h1.savefig(os.path.join(results_dir, - f"convergence_H1_{dim}_nu{nu}.png"), dpi=150) - plt.close(fig_l2); plt.close(fig_h1) - - -# ============================ Main ========================================= -if __name__ == "__main__": - L = 1.0 - E = 1e6 - nx_values = [10, 20, 50, 80, 100, 120, 150, 200, 250] - - element_quad = _QuadElement() - element_tri = _TriElement() - - specs = [ - {"elem": element_quad, "label": "Q1 quad", "marker": "o", "color": "C0"}, - {"elem": element_tri, "label": "P1 tri", "marker": "s", "color": "C1"}, - ] - - DIM = "2d" - for nu in [0.0, 0.3, 0.49, 0.4999]: - print(f"\n nu = {nu} [{DIM}]") - simulation_ponctuelle(element_quad, L, E, nu, nx=20, ny=20, dim=DIM) - simulation_ponctuelle(element_tri, L, E, nu, nx=20, ny=20, dim=DIM) - convergence_study(specs, L, E, nu, nx_values, dim=DIM) - - DIM = "3d" - for nu in [0.0, 0.3, 0.49, 0.4999]: - print(f"\n nu = {nu} [{DIM}]") - simulation_ponctuelle(element_quad, L, E, nu, nx=20, ny=20, dim=DIM) - simulation_ponctuelle(element_tri, L, E, nu, nx=20, ny=20, dim=DIM) - convergence_study(specs, L, E, nu, nx_values, dim=DIM) \ No newline at end of file diff --git a/examples/Freefem/MMS/2D/incompressible.py b/examples/Freefem/MMS/2D/incompressible.py new file mode 100644 index 00000000..280883f9 --- /dev/null +++ b/examples/Freefem/MMS/2D/incompressible.py @@ -0,0 +1,102 @@ +""" +Incompressible 2D MMS on [0,L]^2 with linear-elasticity constitutive law: + + u_ex(x, y) = ( sin(pi x / L) cos(pi y / L), + -cos(pi x / L) sin(pi y / L) ) + + div u = 0 everywhere — the lambda * tr(eps) term in sigma vanishes + pointwise, so the discretization is exercised in the lambda → infinity + (nu → 0.5) limit. Useful for probing volumetric locking. + + BCs: minimal Dirichlet (corner anchor + rotation-killer) so the system is + constrained only by the Neumann tractions on the four edges — that's the + point of the test. +""" + +import numpy as np + +from manufactured_solution import MMSCase2D, lame + + +class Incompressible(MMSCase2D): + name = "incompressible" + plot_label = (r"$u_x = \sin(\pi x/L)\cos(\pi y/L),\ " + r"u_y = -\cos(\pi x/L)\sin(\pi y/L)$") + + def u_ex(self, x, y, L): + k = np.pi / L + return ( np.sin(k * x) * np.cos(k * y), + -np.cos(k * x) * np.sin(k * y)) + + def grad_u_ex(self, x, y, L): + k = np.pi / L + dux_dx = k * np.cos(k * x) * np.cos(k * y) + dux_dy = -k * np.sin(k * x) * np.sin(k * y) + duy_dx = k * np.sin(k * x) * np.sin(k * y) + duy_dy = -k * np.cos(k * x) * np.cos(k * y) + return np.array([[dux_dx, dux_dy], + [duy_dx, duy_dy]]) + + def source(self, x, y, E, nu, L, dim): + lam, mu = lame(E, nu, dim) + k = np.pi / L + ux = np.sin(k * x) * np.cos(k * y) + uy = -np.cos(k * x) * np.sin(k * y) + d2ux_dxx = -k**2 * ux + d2ux_dyy = -k**2 * ux + d2ux_dxy = -k**2 * np.cos(k * x) * np.sin(k * y) + d2uy_dxx = -k**2 * uy + d2uy_dyy = -k**2 * uy + d2uy_dxy = k**2 * np.sin(k * x) * np.cos(k * y) + fx = -((lam + 2*mu) * d2ux_dxx + lam * d2uy_dxy + + mu * (d2ux_dyy + d2uy_dxy)) + fy = -(mu * (d2ux_dxy + d2uy_dxx) + lam * d2ux_dxy + + (lam + 2*mu) * d2uy_dyy) + return (fx, fy) + + def apply_bcs(self, Beam, nodes_2d, L, dim): + # Minimal Dirichlet: pin (0,0) in both directions to kill translation, + # pin uy at (L,0) to kill rotation. Everything else is Neumann. + tmpl = "Vec3d" if dim == "3d" else "Vec2d" + xy = nodes_2d[:, :2] + eps = 1e-10 + idx_corner = [k for k, (xk, yk) in enumerate(xy) + if xk < eps and yk < eps] + idx_rot = [k for k, (xk, yk) in enumerate(xy) + if xk > L - eps and yk < eps] + + if dim == "2d": + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_corner", template=tmpl, + indices=" ".join(map(str, idx_corner)), + fixedDirections="1 1") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_rot", template=tmpl, + indices=" ".join(map(str, idx_rot)), + fixedDirections="0 1") + else: + all_idx = " ".join(map(str, range(len(nodes_2d)))) + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_corner", template=tmpl, + indices=" ".join(map(str, idx_corner)), + fixedDirections="1 1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_rot", template=tmpl, + indices=" ".join(map(str, idx_rot)), + fixedDirections="0 1 0") + Beam.addObject("PartialFixedProjectiveConstraint", + name="fix_z", template=tmpl, + indices=all_idx, + fixedDirections="0 0 1") + + +from beam import case_scene, element_quad, element_tri + +mms = Incompressible() +createScene = case_scene(mms, element_quad) + + +if __name__ == "__main__": + from beam import run_reference_scene + + run_reference_scene(element_quad, mms) diff --git a/examples/Freefem/MMS/2D/params.json b/examples/Freefem/MMS/2D/params.json index 66b001d2..f98c79a1 100644 --- a/examples/Freefem/MMS/2D/params.json +++ b/examples/Freefem/MMS/2D/params.json @@ -7,11 +7,12 @@ "dim": "2d" }, "convergence": { - "nu_values": [0.0, 0.3, 0.49], + "nu_values": [0.0, 0.3, 0.49, 0.4999], "dim_values": ["2d", "3d"], "nx_values": { - "cubic": [10, 20, 30, 40, 50], - "trigonometric": [10, 20, 30, 40, 50] + "cubic": [10, 20, 30, 40, 50], + "trigonometric": [10, 20, 30, 40, 50], + "incompressible": [10, 20, 30, 40, 50] } } } diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py index ad8ad0a5..f877c2c1 100644 --- a/examples/Freefem/MMS/2D/run_convergence.py +++ b/examples/Freefem/MMS/2D/run_convergence.py @@ -8,8 +8,9 @@ import numpy as np import matplotlib.pyplot as plt -from cubic import mms as cubic_mms -from trigonometric import mms as trig_mms +from cubic import mms as cubic_mms +from trigonometric import mms as trig_mms +from incompressible import mms as incomp_mms from beam import ( RESULTS_DIR, @@ -137,7 +138,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): + for mms in (cubic_mms, trig_mms, incomp_mms): nx_vals = conv["nx_values"][mms.name] print(f"\n══ {mms.name} ══") for DIM in conv["dim_values"]: From 0879c43261b4e9bfa014f3c75388dd45a9025f69 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 26 May 2026 15:33:13 +0200 Subject: [PATCH 4/4] Parameterize quadrature rule --- examples/Freefem/MMS/2D/beam.py | 22 ++++++++++++++++--- examples/Freefem/MMS/2D/cubic.py | 10 +++++---- examples/Freefem/MMS/2D/incompressible.py | 10 +++++---- .../Freefem/MMS/2D/manufactured_solution.py | 5 +++++ examples/Freefem/MMS/2D/trigonometric.py | 10 +++++---- 5 files changed, 42 insertions(+), 15 deletions(-) diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py index 4fe58f82..641f593b 100644 --- a/examples/Freefem/MMS/2D/beam.py +++ b/examples/Freefem/MMS/2D/beam.py @@ -86,7 +86,7 @@ def compute_nodal_forces(cls, nodes_2d, conn, mms, L, E, nu, nx, ny, dim): F = assemble_nodal_forces_2d( lambda x, y: mms.source(x, y, E, nu, L, dim), - xy, conn, cls.ELEMENT_RULE) + xy, conn, cls._source_rule(mms)) bottom, top, left, right = _boundary_edges(nx, ny) sides = [(bottom, 0.0, -1.0), @@ -117,7 +117,15 @@ def compute_h1(cls, sol, mms, L): class _QuadElement(_ElementBase): LABEL = "Q1 quad" - ELEMENT_RULE = staticmethod(quad_q1_rule(2)) + 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): @@ -142,7 +150,15 @@ def to_triangles(conn): class _TriElement(_ElementBase): LABEL = "P1 tri" - ELEMENT_RULE = staticmethod(tri_p1_rule(3)) + 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): diff --git a/examples/Freefem/MMS/2D/cubic.py b/examples/Freefem/MMS/2D/cubic.py index 7c94168b..e89e510f 100644 --- a/examples/Freefem/MMS/2D/cubic.py +++ b/examples/Freefem/MMS/2D/cubic.py @@ -11,12 +11,18 @@ import numpy as np from manufactured_solution import MMSCase2D, lame +from beam import (case_scene, run_reference_scene, + element_quad, element_tri, + quad_q1_rule, tri_p1_rule) class Cubic(MMSCase2D): name = "cubic" plot_label = r"$u_x = x(L-x)(x+y),\ u_y = y(L-y)(y-x)$" + source_quadrature_quad = staticmethod(quad_q1_rule(2)) + source_quadrature_tri = staticmethod(tri_p1_rule(3)) + def u_ex(self, x, y, L): return (x * (L - x) * (x + y), y * (L - y) * (y - x)) @@ -75,13 +81,9 @@ def source(self, x, y, E, nu, L, dim): return (fx, fy) -from beam import case_scene, element_quad, element_tri - mms = Cubic() createScene = case_scene(mms, element_quad) if __name__ == "__main__": - from beam import run_reference_scene - run_reference_scene(element_quad, mms) diff --git a/examples/Freefem/MMS/2D/incompressible.py b/examples/Freefem/MMS/2D/incompressible.py index 280883f9..56e1d3ab 100644 --- a/examples/Freefem/MMS/2D/incompressible.py +++ b/examples/Freefem/MMS/2D/incompressible.py @@ -16,6 +16,9 @@ import numpy as np from manufactured_solution import MMSCase2D, lame +from beam import (case_scene, run_reference_scene, + element_quad, element_tri, + quad_q1_rule, tri_p1_rule) class Incompressible(MMSCase2D): @@ -23,6 +26,9 @@ class Incompressible(MMSCase2D): plot_label = (r"$u_x = \sin(\pi x/L)\cos(\pi y/L),\ " r"u_y = -\cos(\pi x/L)\sin(\pi y/L)$") + source_quadrature_quad = staticmethod(quad_q1_rule(2)) + source_quadrature_tri = staticmethod(tri_p1_rule(3)) + def u_ex(self, x, y, L): k = np.pi / L return ( np.sin(k * x) * np.cos(k * y), @@ -90,13 +96,9 @@ def apply_bcs(self, Beam, nodes_2d, L, dim): fixedDirections="0 0 1") -from beam import case_scene, element_quad, element_tri - mms = Incompressible() createScene = case_scene(mms, element_quad) if __name__ == "__main__": - from beam import run_reference_scene - run_reference_scene(element_quad, mms) diff --git a/examples/Freefem/MMS/2D/manufactured_solution.py b/examples/Freefem/MMS/2D/manufactured_solution.py index 573e9bce..a03db6fe 100644 --- a/examples/Freefem/MMS/2D/manufactured_solution.py +++ b/examples/Freefem/MMS/2D/manufactured_solution.py @@ -22,6 +22,11 @@ class MMSCase2D(ABC): name = None # case identifier (must match the params.json key) plot_label = None # LaTeX label for the exact solution + # 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)) + source_quadrature_tri = None # element rule for P1 triangles (e.g. tri_p1_rule(3)) + @abstractmethod def u_ex(self, x, y, L): """Exact solution: returns (ux, uy).""" diff --git a/examples/Freefem/MMS/2D/trigonometric.py b/examples/Freefem/MMS/2D/trigonometric.py index e9527a33..a79400d8 100644 --- a/examples/Freefem/MMS/2D/trigonometric.py +++ b/examples/Freefem/MMS/2D/trigonometric.py @@ -11,6 +11,9 @@ import numpy as np from manufactured_solution import MMSCase2D, lame +from beam import (case_scene, run_reference_scene, + element_quad, element_tri, + quad_q1_rule, tri_p1_rule) class Trigonometric(MMSCase2D): @@ -18,6 +21,9 @@ class Trigonometric(MMSCase2D): plot_label = (r"$u_x = \sin(\pi x/L)\cos(\pi y/L),\ " r"u_y = \cos(\pi x/L)\sin(\pi y/L)$") + source_quadrature_quad = staticmethod(quad_q1_rule(2)) + source_quadrature_tri = staticmethod(tri_p1_rule(3)) + def u_ex(self, x, y, L): k = np.pi / L return (np.sin(k * x) * np.cos(k * y), @@ -81,13 +87,9 @@ def source(self, x, y, E, nu, L, dim): return (fx, fy) -from beam import case_scene, element_quad, element_tri - mms = Trigonometric() createScene = case_scene(mms, element_quad) if __name__ == "__main__": - from beam import run_reference_scene - run_reference_scene(element_quad, mms)