{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ec11fe4",
   "metadata": {},
   "source": [
    "# GaLS MHM: write displacement–pressure equations with UFL\n",
    "\n",
    "This tutorial builds the complete method from local and global equations. It uses the unit square, a smooth displacement, shear modulus $\\mu=1$ and Lamé modulus $\\lambda=4999$. The classical control is conforming Taylor–Hood on independent global meshes. See [Gomes, Pereira and Valentin (2024)](https://arxiv.org/abs/2403.16890v1) for the stabilized formulation and its hypotheses.\n",
    "\n",
    "Run in the locked `introduction` Pixi environment with DOLFINx, UFL and Basix. The source notebook and its companion examples are separate from the installed PyMHM package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "companion-bootstrap",
   "metadata": {
    "tags": [
     "pymhm-companion-bootstrap"
    ]
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import os\n",
    "import sys\n",
    "from pymhm.io.workspace import workspace_from_archive\n",
    "\n",
    "# Download checksum-verified support files; extraction does not execute them.\n",
    "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/4815ac8e08194cc102af53a79bb5e19a8c536656aed1634e1faa3417ed816b67/gals_mhm_workflow-companion.zip\"\n",
    "COMPANION_SHA256 = \"4815ac8e08194cc102af53a79bb5e19a8c536656aed1634e1faa3417ed816b67\"\n",
    "WORKSPACE = Path(\n",
    "    os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd() / \".pymhm-companions\" / COMPANION_SHA256)\n",
    ")\n",
    "ROOT = workspace_from_archive(COMPANION_URL, sha256=COMPANION_SHA256, directory=WORKSPACE)\n",
    "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n",
    "if str(ROOT) not in sys.path:\n",
    "    sys.path.insert(0, str(ROOT))\n",
    "\n",
    "# Prepare only this notebook's declared inputs, without executing its equations.\n",
    "from scripts.notebook_reproduction import notebook_workspace\n",
    "\n",
    "ROOT = notebook_workspace(\"elasticity/gals_mhm_workflow.ipynb\", directory=ROOT)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3ccfbd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "ROOT = Path(os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd())).resolve()\n",
    "if not (ROOT / \"examples\").is_dir():\n",
    "    ROOT = next(parent for parent in ROOT.parents if (parent / \"examples\").is_dir())\n",
    "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n",
    "if str(ROOT) not in sys.path:\n",
    "    sys.path.insert(0, str(ROOT))\n",
    "print(\"Workspace:\", ROOT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77554205",
   "metadata": {},
   "outputs": [],
   "source": [
    "from types import SimpleNamespace\n",
    "\n",
    "import basix.ufl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import ufl\n",
    "from threadpoolctl import threadpool_limits\n",
    "\n",
    "from pymhm import (\n",
    "    Equation,\n",
    "    FaceSpace,\n",
    "    MeshHierarchy,\n",
    "    MultiscaleProblem,\n",
    "    SkeletonSpace,\n",
    "    TriangleMesh,\n",
    "    assemble,\n",
    "    bind_interface,\n",
    "    bind_problem,\n",
    "    columns,\n",
    "    solve,\n",
    ")\n",
    "from pymhm.core.equations import compile_form\n",
    "from pymhm.fem.scalar.triangle import nodal_space\n",
    "from pymhm.fem.vector.elasticity import rigid_modes\n",
    "from pymhm.postprocessing.solutions import ElasticitySolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46fc4263",
   "metadata": {},
   "source": [
    "## 1. Define physical data before choosing the method\n",
    "\n",
    "The sign convention is\n",
    "\n",
    "$$\n",
    "p=-\\lambda\\nabla\\cdot u,\\qquad\n",
    "\\sigma=2\\mu\\varepsilon(u)-pI,\\qquad -\\nabla\\cdot\\sigma=f.\n",
    "$$\n",
    "\n",
    "The solenoidal trigonometric displacement has a small compressible correction proportional to $\\beta=\\mu/(\\lambda+\\mu)$. UFL differentiates the exact stress to derive the source. Independent NumPy expressions evaluate the exact fields for errors and figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de079863",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu, lam, alpha = 1.0, 4999.0, 0.25\n",
    "\n",
    "\n",
    "def exact_ufl(domain):\n",
    "    \"\"\"Differentiate the exact Cauchy stress symbolically, independently of assembly.\"\"\"\n",
    "    x = ufl.SpatialCoordinate(domain)\n",
    "    beta = mu / (lam + mu)\n",
    "    ue = ufl.as_vector(\n",
    "        (\n",
    "            (ufl.cos(2 * np.pi * x[0]) - 1) * ufl.sin(2 * np.pi * x[1]),\n",
    "            (1 - ufl.cos(2 * np.pi * x[1])) * ufl.sin(2 * np.pi * x[0]),\n",
    "        )\n",
    "    )\n",
    "    ue += beta * ufl.sin(np.pi * x[0]) * ufl.sin(np.pi * x[1]) * ufl.as_vector((1, 1))\n",
    "    pe = -lam * beta * np.pi * ufl.sin(np.pi * (x[0] + x[1]))\n",
    "    exact_sigma = 2 * mu * ufl.sym(ufl.grad(ue)) - pe * ufl.Identity(2)\n",
    "    force = -ufl.div(exact_sigma)\n",
    "    return ue, pe, force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1600dd55",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_u(points):\n",
    "    \"\"\"Evaluate the boundary-vanishing displacement, including finite compressibility.\"\"\"\n",
    "    x, y = points.T\n",
    "    beta = mu / (lam + mu)\n",
    "    value = np.column_stack(\n",
    "        (\n",
    "            (np.cos(2 * np.pi * x) - 1) * np.sin(2 * np.pi * y),\n",
    "            (1 - np.cos(2 * np.pi * y)) * np.sin(2 * np.pi * x),\n",
    "        )\n",
    "    )\n",
    "    return value + beta * (np.sin(np.pi * x) * np.sin(np.pi * y))[:, None]\n",
    "\n",
    "\n",
    "def exact_p(points):\n",
    "    \"\"\"Evaluate the zero-mean Herrmann pressure dictated by the displacement.\"\"\"\n",
    "    return -lam * mu / (lam + mu) * np.pi * np.sin(np.pi * points.sum(axis=1))\n",
    "\n",
    "\n",
    "def exact_gradient(points):\n",
    "    \"\"\"Differentiate each displacement component in physical coordinates.\"\"\"\n",
    "    x, y = points.T\n",
    "    sx, sy, cx, cy = (\n",
    "        np.sin(2 * np.pi * x),\n",
    "        np.sin(2 * np.pi * y),\n",
    "        np.cos(2 * np.pi * x),\n",
    "        np.cos(2 * np.pi * y),\n",
    "    )\n",
    "    gradient = (\n",
    "        2\n",
    "        * np.pi\n",
    "        * np.stack(\n",
    "            (np.column_stack((-sx * sy, (cx - 1) * cy)), np.column_stack(((1 - cy) * cx, sx * sy))),\n",
    "            axis=1,\n",
    "        )\n",
    "    )\n",
    "    correction = (\n",
    "        mu\n",
    "        / (lam + mu)\n",
    "        * np.pi\n",
    "        * np.column_stack(\n",
    "            (np.cos(np.pi * x) * np.sin(np.pi * y), np.sin(np.pi * x) * np.cos(np.pi * y))\n",
    "        )\n",
    "    )\n",
    "    return gradient + correction[:, None, :]\n",
    "\n",
    "\n",
    "def exact_stress(points):\n",
    "    \"\"\"Evaluate the physical symmetric Cauchy stress.\"\"\"\n",
    "    gradient = exact_gradient(points)\n",
    "    return mu * (gradient + gradient.swapaxes(1, 2)) - exact_p(points)[:, None, None] * np.eye(2)\n",
    "\n",
    "\n",
    "exact = SimpleNamespace(\n",
    "    displacement=exact_u, pressure=exact_p, gradient=exact_gradient, stress=exact_stress\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e27f713",
   "metadata": {},
   "source": [
    "## 2. Choose macro, local and interface spaces\n",
    "\n",
    "There are $128$ macrotriangles. Each has its own refinement $r=4$, continuous vector $P_1$ displacement and continuous $P_1$ pressure. Each macroface carries two components of discontinuous $P_1$ negative traction. This $(k,r)=(1,4)$ choice resolves the linear trace; equal-order spaces require GaLS.\n",
    "\n",
    "`MeshHierarchy` associates the independent local meshes. `bind_interface` handles geometric incidence signs; the physical signs remain in the equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9603385a",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu, lam, alpha = 1.0, 4999.0, 0.25\n",
    "degree, refinement, n = 1, 4, 8\n",
    "macro = TriangleMesh.unit_square(n)\n",
    "skeleton = SkeletonSpace(macro, tuple(FaceSpace.uniform(1) for _ in macro.faces), 2)\n",
    "hierarchy = MeshHierarchy(\n",
    "    macro, tuple(macro.submesh(i, refinement) for i in range(len(macro.cells)))\n",
    ")\n",
    "interface = bind_interface(skeleton, convention=\"normal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "666f1b31",
   "metadata": {},
   "source": [
    "## 3. Translate every term of the local variational problem\n",
    "\n",
    "Let $R(u,p)=\\nabla\\cdot(2\\mu\\varepsilon(u))-\\nabla p$ and $\\tau_h=\\alpha h^2$. The local volume forms are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_T((u,p),(v,q))\n",
    "&=(2\\mu\\varepsilon(u),\\varepsilon(v))_T\n",
    " -(p,\\nabla\\cdot v)_T-(q,\\nabla\\cdot u)_T\\\\\n",
    "&\\quad-\\lambda^{-1}(p,q)_T\n",
    " -\\sum_{t\\subset T}(\\tau_hR(u,p),R(v,q))_t,\\\\\n",
    "L_T(v,q)&=(f,v)_T+\\sum_{t\\subset T}(\\tau_h f,R(v,q))_t.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Both the negative residual product and the positive source pairing are needed. For constant shear and affine displacement elements, the strong displacement residual vanishes inside each fine triangle. The inverse bound gives $0<\\alpha<1/(2\\mu)$; $\\alpha=1/4$ satisfies it here. It must be recomputed for higher degree or variable shear.\n",
    "\n",
    "The local kernel consists of the two translations and one rigid rotation, with pressure zero. Their physical displacement integrals fix the local lifting. `columns` records these moment forms, while `native.to_native` delegates coefficient ordering to the space adapter. `b` pairs negative traction with displacement; `c` independently records the negative transpose balance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "972cce89",
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_equations(local):\n",
    "    \"\"\"Declare native P1/P1 GaLS forms, rigid moments and negative physical traction.\"\"\"\n",
    "    element = basix.ufl.mixed_element(\n",
    "        [\n",
    "            basix.ufl.element(\n",
    "                \"Lagrange\",\n",
    "                \"triangle\",\n",
    "                1,\n",
    "                shape=(2,),\n",
    "                lagrange_variant=basix.LagrangeVariant.equispaced,\n",
    "            ),\n",
    "            basix.ufl.element(\n",
    "                \"Lagrange\", \"triangle\", 1, lagrange_variant=basix.LagrangeVariant.equispaced\n",
    "            ),\n",
    "        ]\n",
    "    )\n",
    "    native = local.native_space(element)\n",
    "    domain = native.space.mesh\n",
    "    x = ufl.SpatialCoordinate(domain)\n",
    "    u, p = ufl.TrialFunctions(native.space)\n",
    "    v, q = ufl.TestFunctions(native.space)\n",
    "    ue, pe, force = exact_ufl(domain)\n",
    "    dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 8})\n",
    "\n",
    "    def eps(w):\n",
    "        \"\"\"Return the symmetric physical displacement gradient.\"\"\"\n",
    "        return ufl.sym(ufl.grad(w))\n",
    "\n",
    "    def residual(w, r):\n",
    "        \"\"\"Return the signed strong divergence of physical Cauchy stress.\"\"\"\n",
    "        return ufl.div(2 * mu * eps(w)) - ufl.grad(r)\n",
    "\n",
    "    tau = alpha * ufl.CellDiameter(domain) ** 2\n",
    "    a = (\n",
    "        2 * mu * ufl.inner(eps(u), eps(v))\n",
    "        - p * ufl.div(v)\n",
    "        - q * ufl.div(u)\n",
    "        - p * q / lam\n",
    "        - tau * ufl.inner(residual(u, p), residual(v, q))\n",
    "    ) * dx\n",
    "    L = (ufl.inner(force, v) + tau * ufl.inner(force, residual(v, q))) * dx\n",
    "    center = macro.points.mean(axis=0)\n",
    "    rigid = [\n",
    "        ufl.as_vector((1, 0)),\n",
    "        ufl.as_vector((0, 1)),\n",
    "        ufl.as_vector((-(x[1] - center[1]), x[0] - center[0])),\n",
    "    ]\n",
    "    nodes = nodal_space(local.mesh, 1)[1]\n",
    "    kernel = native.to_native(rigid_modes(nodes, center).reshape(2 * len(nodes), 3), component=0)\n",
    "    moments = columns(*(ufl.inner(v, z) * dx for z in rigid))\n",
    "    b = local.trace_pairings(lambda multiplier, ds: ufl.inner(multiplier, v) * ds)\n",
    "    c = local.trace_pairings(lambda test_trace, ds: -ufl.inner(test_trace, v) * ds, axis=\"rows\")\n",
    "    pressure_moment = compile_form(q * dx)\n",
    "    local.field(\"displacement\", native, component=0)\n",
    "    local.field(\"pressure\", native, component=1)\n",
    "    return local.equations(\n",
    "        a=a, L=L, b=b, c=c, kernel=kernel, moments=moments, metadata=(pressure_moment,)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f2578d5",
   "metadata": {},
   "source": [
    "## 4. State the global equation and the physical pressure identity\n",
    "\n",
    "The complete displacement is prescribed weakly at the outer boundary:\n",
    "\n",
    "$$\n",
    "\\sum_T\\langle u_T,\\mu\\rangle_{\\partial T}\n",
    "=\\langle g,\\mu\\rangle_{\\partial\\Omega},\\qquad g=0.\n",
    "$$\n",
    "\n",
    "Three rigid coordinates per macrocell remain in the global problem. With finite $\\lambda$, pressure is determined by compressibility:\n",
    "\n",
    "$$\n",
    "\\int_\\Omega p/\\lambda=-\\int_{\\partial\\Omega}g\\cdot n=0.\n",
    "$$\n",
    "\n",
    "For this constant material, the equivalent normalized row is $\\int_\\Omega p=0$. It is a physical identity, not a freely selected pressure gauge. At $\\lambda=\\infty$, zero net boundary displacement is necessary and the pressure mean becomes a gauge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93413cf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def global_equation(global_context):\n",
    "    \"\"\"Impose homogeneous weak displacement through the global trace functional.\"\"\"\n",
    "    boundary, _ = global_context.boundary_data((0, 0), order=8)\n",
    "    return Equation(0, global_context.trace_load(-boundary))\n",
    "\n",
    "\n",
    "with threadpool_limits(1):\n",
    "    problem = bind_problem(\n",
    "        hierarchy, interface, local_equations, global_equation=global_equation, retained=3\n",
    "    )\n",
    "    system = assemble(problem)\n",
    "    pressure_identity = system.mean_constraint([entry[0] for entry in system.local_metadata], 0.0)\n",
    "    coefficients = solve(system, constraints=[pressure_identity])\n",
    "    ufields = coefficients.field(\"displacement\")\n",
    "    pfields = coefficients.field(\"pressure\")\n",
    "    view = ElasticitySolution(\n",
    "        skeleton,\n",
    "        tuple(f.mesh for f in ufields),\n",
    "        tuple(f.portable_coefficients.reshape(-1, 2) for f in ufields),\n",
    "        tuple(f.portable_coefficients for f in pfields),\n",
    "        coefficients,\n",
    "        1,\n",
    "        1,\n",
    "        lam,\n",
    "        mu,\n",
    "        \"gals\",\n",
    "        (alpha,) * len(ufields),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c874d07d",
   "metadata": {},
   "source": [
    "## 5. Recover physical fields and test the original equations\n",
    "\n",
    "The field declarations carry their executed native maps into the recovered solution. The `ElasticitySolution` view applies the physical constitutive law to those fields; it does not choose or solve a method. We check displacement and pressure equations separately before combining their original row residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99b9d425",
   "metadata": {},
   "outputs": [],
   "source": [
    "from examples.core_elasticity_field_archive import ProductionObservation\n",
    "from examples.minimal_flow_originals import original_diagnostics\n",
    "\n",
    "errors = {\n",
    "    \"displacement_l2\": view.l2_error(exact_u, 12),\n",
    "    \"pressure_l2\": view.pressure_l2_error(exact_p, 12),\n",
    "    \"stress_l2\": view.stress_l2_error(exact_stress, 12),\n",
    "    \"gradient_l2\": view.h1_seminorm_error(exact_gradient, 12),\n",
    "}\n",
    "print(errors)\n",
    "print(\"Reduced residual:\", coefficients.residual)\n",
    "\n",
    "observation = ProductionObservation(\n",
    "    system=system,\n",
    "    applied_boundary=-system.global_load[: system.trace_size],\n",
    "    constraints=[pressure_identity],\n",
    ")\n",
    "blocks = [\n",
    "    {\n",
    "        \"displacement_equilibrium\": len(field.portable_coefficients),\n",
    "        \"pressure_compressibility\": len(pfield.portable_coefficients),\n",
    "    }\n",
    "    for field, pfield in zip(ufields, pfields, strict=True)\n",
    "]\n",
    "physical = original_diagnostics(view, observation, blocks)\n",
    "print(\n",
    "    \"Original uncondensed physical residual:\", physical[\"full_uncondensed_relative_to_physical_rhs\"]\n",
    ")\n",
    "assert physical[\"accepted\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc27017a",
   "metadata": {},
   "source": [
    "## 6. Assemble an independent fine classical reference\n",
    "\n",
    "The following function writes the unstabilized conforming Taylor–Hood $P_2/P_1$ formulation. Its displacement boundary values are imposed strongly on global meshes with $32$ and $64$ squares per axis. The same material, source and pressure identity are used. Its own analytical errors must decrease before the fine field is used as a baseline. These global meshes are not passed to MHM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b14c800",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfinx import fem\n",
    "\n",
    "from pymhm.backends.spaces import bind_space, coefficient_map\n",
    "from pymhm.postprocessing.fields import DiscreteField, FieldDefinition\n",
    "\n",
    "\n",
    "def conforming_reference(resolution):\n",
    "    \"\"\"Solve the notebook's unstabilized Taylor–Hood form on a fine global mesh.\"\"\"\n",
    "    mesh = TriangleMesh.unit_square(resolution)\n",
    "    element = basix.ufl.mixed_element(\n",
    "        [\n",
    "            basix.ufl.element(\n",
    "                \"Lagrange\",\n",
    "                \"triangle\",\n",
    "                2,\n",
    "                shape=(2,),\n",
    "                lagrange_variant=basix.LagrangeVariant.equispaced,\n",
    "            ),\n",
    "            basix.ufl.element(\n",
    "                \"Lagrange\", \"triangle\", 1, lagrange_variant=basix.LagrangeVariant.equispaced\n",
    "            ),\n",
    "        ]\n",
    "    )\n",
    "    native = bind_space(mesh, element)\n",
    "    try:\n",
    "        u, p = ufl.TrialFunctions(native.space)\n",
    "        v, q = ufl.TestFunctions(native.space)\n",
    "        domain = native.space.mesh\n",
    "        ue, pe, force = exact_ufl(domain)\n",
    "        dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 12})\n",
    "        a = (\n",
    "            2 * mu * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(v)))\n",
    "            - p * ufl.div(v)\n",
    "            - q * ufl.div(u)\n",
    "            - p * q / lam\n",
    "        ) * dx\n",
    "        L = ufl.inner(force, v) * dx\n",
    "        _, nodes = nodal_space(mesh, 2)\n",
    "        boundary_nodes = np.flatnonzero(np.any(np.isclose(nodes, 0) | np.isclose(nodes, 1), axis=1))\n",
    "        portable = (2 * boundary_nodes[:, None] + np.arange(2)).ravel()\n",
    "        native_ids = coefficient_map(native, component=0)[portable]\n",
    "        fixed = dict(zip(native_ids, exact_u(nodes[boundary_nodes]).ravel(), strict=True))\n",
    "        reference_problem = MultiscaleProblem.from_global(\n",
    "            Equation(a, L),\n",
    "            native.size,\n",
    "            fixed=fixed,\n",
    "        )\n",
    "        # Finite lambda determines this zero integral by the displacement data.\n",
    "        pressure_integral = compile_form(q * dx)\n",
    "        result = solve(assemble(reference_problem), constraints=[(pressure_integral, 0.0)])\n",
    "        numerical = fem.Function(native.space)\n",
    "        numerical.x.array[:] = result.trace\n",
    "        uh, ph = ufl.split(numerical)\n",
    "        sigma_h = 2 * mu * ufl.sym(ufl.grad(uh)) - ph * ufl.Identity(2)\n",
    "        sigma_e = 2 * mu * ufl.sym(ufl.grad(ue)) - pe * ufl.Identity(2)\n",
    "        differences = {\n",
    "            \"displacement_l2\": uh - ue,\n",
    "            \"pressure_l2\": ph - pe,\n",
    "            \"stress_l2\": sigma_h - sigma_e,\n",
    "        }\n",
    "        errors = {\n",
    "            name: float(np.sqrt(fem.assemble_scalar(fem.form(ufl.inner(delta, delta) * dx))))\n",
    "            for name, delta in differences.items()\n",
    "        }\n",
    "        fields = {\n",
    "            name: DiscreteField(\n",
    "                FieldDefinition(name, descriptor=native.descriptor(component=i)), result.trace\n",
    "            )\n",
    "            for i, name in enumerate((\"displacement\", \"pressure\"))\n",
    "        }\n",
    "        return fields, {\n",
    "            \"resolution\": resolution,\n",
    "            \"triangles\": len(mesh.cells),\n",
    "            \"unknowns\": native.size,\n",
    "            \"errors\": errors,\n",
    "            \"residual\": result.residual,\n",
    "        }\n",
    "    finally:\n",
    "        native.close()\n",
    "\n",
    "\n",
    "with threadpool_limits(1):\n",
    "    reference32, control32 = conforming_reference(32)\n",
    "    reference64, control64 = conforming_reference(64)\n",
    "print(control32)\n",
    "print(control64)\n",
    "assert all(\n",
    "    control64[\"errors\"][name] < control32[\"errors\"][name]\n",
    "    for name in (\"displacement_l2\", \"stress_l2\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3863439a",
   "metadata": {},
   "source": [
    "## 7. Plot the physical fields on their actual meshes\n",
    "\n",
    "The display helper preserves independent one-sided MHM samples and draws the macro mesh on every panel. Each panel has an independent color scale. The coarse MHM field and fine conforming field solve the same operator in different spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927d2206",
   "metadata": {},
   "outputs": [],
   "source": [
    "from examples.introduction.vector import plot_field_panels\n",
    "from examples.plot_elasticity import field_arrays\n",
    "\n",
    "samples = field_arrays(view, exact)\n",
    "x, triangles = samples[\"points\"], samples[\"cells\"]\n",
    "baseline_u = reference64[\"displacement\"].evaluate(x)\n",
    "panels = {\n",
    "    \"Exact displacement\\nmagnitude\": (\n",
    "        x,\n",
    "        triangles,\n",
    "        np.linalg.norm(samples[\"exact_displacement\"], axis=1),\n",
    "    ),\n",
    "    \"GaLS displacement\\nmagnitude\": (x, triangles, np.linalg.norm(samples[\"displacement\"], axis=1)),\n",
    "    \"Fine Taylor–Hood\\ndisplacement magnitude\": (x, triangles, np.linalg.norm(baseline_u, axis=1)),\n",
    "    \"Exact Herrmann\\npressure\": (x, triangles, samples[\"exact_pressure\"]),\n",
    "    \"GaLS Herrmann\\npressure\": (x, triangles, samples[\"pressure\"]),\n",
    "    \"GaLS displacement\\nerror\": (\n",
    "        x,\n",
    "        triangles,\n",
    "        np.linalg.norm(samples[\"displacement\"] - samples[\"exact_displacement\"], axis=1),\n",
    "    ),\n",
    "}\n",
    "figure = plot_field_panels(macro, panels, figsize=(13, 8))\n",
    "output = ROOT / \"build\" / \"tutorial-method-workflows\"\n",
    "output.mkdir(parents=True, exist_ok=True)\n",
    "figure.savefig(output / \"gals-elasticity-fields.png\", dpi=180, bbox_inches=\"tight\")\n",
    "figure.savefig(output / \"gals-elasticity-fields.svg\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1790bc80",
   "metadata": {},
   "source": [
    "## 8. Refine the skeletal and local discretizations\n",
    "\n",
    "The smooth convergence campaign uses this same trigonometric data and $\\lambda/\\mu=4999$. It holds the $4\\times4$ macro grid fixed, subdivides every trace into $s=1,2,3,4,6,8$ segments and chooses local refinement $r=4s$, so the fine-element and trace-segment sizes decrease together. The plotted illustration above uses an independent $8\\times8$ macro grid with one trace segment per face. The campaign's rows are attributed archived acquisitions, not recomputed by reading them.\n",
    "\n",
    "For this smooth admissible $P_1/P_1$ family, the reference order is one for the\n",
    "broken gradient. Theorem 5.2 of\n",
    "[Gomes, Pereira and Valentin (2024)](https://arxiv.org/abs/2403.16890v1)\n",
    "provides the energy estimate. Displacement $L^2$ order two additionally uses\n",
    "the local and global dual smoothing assumptions of Remark 5.5; it does not\n",
    "follow from that energy theorem alone. The smooth square-domain study keeps\n",
    "the constant material and compatible approximation family fixed and measures\n",
    "both norms separately. Equation (5.8) in Theorem 5.2 also includes the pressure $L^2$ and stress $L^2$ errors in its combined order-one upper bound. An upper bound does not require every component to have exactly that observed exponent: the measured pressure orders are stronger and are reported separately.\n",
    "Changing only the macro grid or only the trace with a fixed local resolution\n",
    "defines another study.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64b63770",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "from examples.tutorial_convergence import refinement_series,plot_method_series,asymptotic_summary\n",
    "record_path=ROOT/'examples/results/elasticity.json'\n",
    "record=json.loads(record_path.read_text())\n",
    "rows=[row for row in record['refinement'] if row['method']=='gals-p1']\n",
    "series=refinement_series(\n",
    "    record_path,rows,'skeleton_size',\n",
    "    {'displacement L2':'displacement_l2','pressure L2':'pressure_l2','broken gradient L2':'gradient_l2'},\n",
    "    {'displacement L2':2.0,'broken gradient L2':1.0},method='gals-elasticity',\n",
    "    spaces='P1/P1, segmented P1 traction; local refinement four times segment count',\n",
    "    rate_provenance='Gomes, Pereira and Valentin (2024), theorem 5.2: conditional smooth energy order one; displacement L2 order two additionally requires Remark 5.5 dual smoothing; pressure orders reported separately',\n",
    "    refinement='trace/fine size on a fixed macro mesh',root=ROOT,\n",
    ")\n",
    "print(asymptotic_summary(series))\n",
    "figure=plot_method_series(series)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d2bfcf8",
   "metadata": {},
   "source": [
    "The trigonometric data here are an analytical verification with explicitly derived source; they are not a claim of matching every table or local solver in the paper.\n",
    "\n",
    "## References\n",
    "\n",
    "- [Gomes, Pereira and Valentin (2024), *A low-order locking-free multiscale finite element method for isotropic elasticity*, preprint v1](https://arxiv.org/abs/2403.16890v1).\n",
    "- [Harder, Madureira and Valentin (2016), *A hybrid-mixed method for elasticity*](https://doi.org/10.1051/m2an/2015046).\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "pymhm": {
   "problem": "elasticity",
   "title": "GaLS MHM through explicit UFL local and global equations",
   "workflow": "mesh-first explicit local/global equations"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}