{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0e935eb1",
   "metadata": {},
   "source": [
    "# Independent skeletal refinement for unfitted MHM\n",
    "\n",
    "An unfitted trace partition can cut across the local boundary elements. It is declared independently rather than forcing its breaks to coincide with the fine mesh. We first show this geometric capability on smooth identity diffusion with an exact pressure, then distinguish its qualified convergence study from material-interface cases.\n",
    "\n",
    "The flux-approximation theorem and trace construction follow [Chaumont-Frelet, Paredes and Valentin (2026)](https://doi.org/10.1016/j.camwa.2026.01.016). This smooth example is not a reproduction of their heterogeneous applications."
   ]
  },
  {
   "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/4605b8257729af9879b5e9294568a59c3b8253d50adbb3bb7e7fe7a77ff5d149/unfitted_trace_workflow-companion.zip\"\n",
    "COMPANION_SHA256 = \"4605b8257729af9879b5e9294568a59c3b8253d50adbb3bb7e7fe7a77ff5d149\"\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(\"darcy/unfitted_trace_workflow.ipynb\", directory=ROOT)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09c57ecf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ufl\n",
    "from pymhm import (\n",
    "    Equation, FaceSpace, LocalContext, MeshHierarchy, SkeletonSpace, TriangleMesh,\n",
    "    assemble, bind_interface, bind_problem, columns, solve,\n",
    ")\n",
    "from threadpoolctl import threadpool_limits\n",
    "\n",
    "\n",
    "def exact_pressure(points: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Analytical pressure; it vanishes on every exterior edge.\"\"\"\n",
    "    return np.sin(2 * np.pi * points[:, 0]) * np.sin(2 * np.pi * points[:, 1])\n",
    "\n",
    "\n",
    "def exact_gradient(points: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Independently differentiated physical pressure gradient.\"\"\"\n",
    "    x, y = 2 * np.pi * points.T\n",
    "    return 2 * np.pi * np.column_stack((np.cos(x) * np.sin(y), np.sin(x) * np.cos(y)))\n",
    "\n",
    "\n",
    "def source(points: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Negative Laplacian computed independently of the discrete operator.\"\"\"\n",
    "    return 8 * np.pi**2 * exact_pressure(points)\n",
    "from pymhm.fem.scalar.triangle import scalar_operators, trace_coupling\n",
    "from pymhm.postprocessing.nodal import nodal_field\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06426526",
   "metadata": {},
   "source": [
    "## 1. Define macro, local and interface spaces\n",
    "\n",
    "The hierarchy owns mesh association. The interface binding owns geometric incidence and numbering; variational signs are declared separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c31ebe04",
   "metadata": {},
   "outputs": [],
   "source": [
    "macro = TriangleMesh.unit_square(2)\n",
    "local_meshes = tuple(macro.submesh(cell, 4) for cell in range(len(macro.cells)))\n",
    "hierarchy = MeshHierarchy(macro, local_meshes)\n",
    "skeleton = SkeletonSpace(macro, tuple(FaceSpace.uniform(1, subdivisions=3) for _ in macro.faces))\n",
    "interface = bind_interface(skeleton, convention=\"normal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8edff754",
   "metadata": {},
   "source": [
    "## 2. Write the local primal form\n",
    "\n",
    "$$\n",
    "(\\nabla p_T,\\nabla v)_T+\\langle\\lambda_T,v\\rangle_{\\partial T}=(f,v)_T.\n",
    "$$\n",
    "\n",
    "Each macroedge has three P1 trace segments; each local edge has four fine subdivisions with P3 volume pressure. Their partitions differ. The portable Basix owner `scalar_operators` integrates the declared volume form; `trace_coupling` integrates on the common intersections and applies canonical incidence. Its global coordinates are declared explicitly. The UFL volume integrand would be `ufl.inner(ufl.grad(p), ufl.grad(v)) * dx`; native `trace_pairings` currently requires boundary facets aligned with trace breaks, so it is not used for this independent partition. The constant kernel and physical volume moment remain explicit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df2b4ea0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_equations(local: LocalContext):\n",
    "    \"\"\"Declare volume diffusion, physical mean and signed scalar trace forms.\"\"\"\n",
    "    fine = local.mesh\n",
    "    A, mass, force = scalar_operators(fine, 3, diffusion=1.0, source=source, order=12)\n",
    "    B = trace_coupling(macro, local.cell, fine, skeleton, 3)\n",
    "    constant = np.ones((A.shape[0], 1))\n",
    "    field = nodal_field(\"pressure\", fine, 3)\n",
    "    local.field(\"pressure\", evaluator=field.evaluator,\n",
    "                gradient_evaluator=field.gradient_evaluator, basis_id=field.basis_id)\n",
    "    return local.equations(a=A, L=force, b=B, c=-B.T,\n",
    "                          kernel=constant, moments=mass @ constant,\n",
    "                          coordinates=\"global\", metadata=(fine,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be5eaf8f",
   "metadata": {},
   "source": [
    "## 3. State global pressure continuity\n",
    "\n",
    "$$\n",
    "-\\sum_T\\langle p_T,\\mu_T\\rangle_{\\partial T}=-\\langle p_D,\\mu\\rangle_{\\partial\\Omega}.\n",
    "$$\n",
    "\n",
    "Retain one physical pressure mean per macrocell. A material-interface study additionally specifies region markers, each region's coefficient, physical continuity and certified bounds. A trace partition by itself does not supply those data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc1d7be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def global_equation(global_context):\n",
    "    \"\"\"Supply the weak Dirichlet moments in the bound interface layout.\"\"\"\n",
    "    boundary, fixed = global_context.boundary_data(0.0, order=8)\n",
    "    assert not fixed\n",
    "    return Equation(0, global_context.trace_load(-boundary))\n",
    "\n",
    "\n",
    "problem = bind_problem(hierarchy, interface, local_equations,\n",
    "                       global_equation=global_equation, retained=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70d5bd60",
   "metadata": {},
   "source": [
    "## 4. Assemble and solve\n",
    "\n",
    "Assembly solves the independent local equations and sums their interface contributions. Solving reconstructs named local fields in their executed coefficient bases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21963788",
   "metadata": {},
   "outputs": [],
   "source": [
    "with threadpool_limits(1):\n",
    "    system = assemble(problem)\n",
    "    coefficients = solve(system)\n",
    "pressure_fields = coefficients.field(\"pressure\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "debcc177",
   "metadata": {},
   "source": [
    "## 5. Integrate physical errors\n",
    "\n",
    "Pressure and its broken gradient are integrated separately with positive order-12 quadrature. No algebraic residual is used as a substitute for these errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d24f65ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymhm.fem.scalar.operators import triangle_quadrature\n",
    "bary, weights = triangle_quadrature(12)\n",
    "pressure_error, gradient_error = 0.0, 0.0\n",
    "for fine, field in zip(local_meshes, pressure_fields, strict=True):\n",
    "    points = np.einsum(\"qi,tia->tqa\", bary, fine.points[fine.cells])\n",
    "    flat = points.reshape(-1, 2)\n",
    "    ep = (field.evaluate(flat) - exact_pressure(flat)).reshape(points.shape[:2])\n",
    "    eg = (field.gradient(flat) - exact_gradient(flat)).reshape(*points.shape[:2], 2)\n",
    "    pressure_error += np.einsum(\"t,q,tq->\", fine.areas, weights, ep**2)\n",
    "    gradient_error += np.einsum(\"t,q,tq->\", fine.areas, weights, np.sum(eg**2, axis=-1))\n",
    "print({\"original_equation_residual\": coefficients.raw_residual,\n",
    "       \"pressure_L2\": float(np.sqrt(pressure_error)),\n",
    "       \"broken_gradient_L2\": float(np.sqrt(gradient_error))})\n",
    "assert coefficients.raw_residual < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a538df",
   "metadata": {},
   "source": [
    "## 6. Plot analytical, numerical and error fields\n",
    "\n",
    "Independent local panels retain one-sided values. Every panel shows the actual macro partition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a4f604a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(8.5, 7), layout=\"constrained\")\n",
    "all_exact, all_values, panels = [], [], []\n",
    "for fine, field in zip(local_meshes, pressure_fields, strict=True):\n",
    "    points = fine.points[fine.cells].mean(axis=1)\n",
    "    exact, values = exact_pressure(points), field.evaluate(points)\n",
    "    flux_magnitude = np.linalg.norm(-field.gradient(points), axis=1)\n",
    "    panels.append((fine, exact, values, flux_magnitude))\n",
    "    all_exact.extend(exact)\n",
    "    all_values.extend(values)\n",
    "common = max(np.max(np.abs(all_exact)), np.max(np.abs(all_values)))\n",
    "error_max = max(np.max(abs(values - exact)) for _, exact, values, _ in panels)\n",
    "flux_max = max(np.max(flux) for _, _, _, flux in panels)\n",
    "titles = (\"Analytical pressure\", \"Hybrid pressure\", \"Pressure error\", \"Hybrid Darcy flux magnitude\")\n",
    "for index, (ax, title) in enumerate(zip(axes.flat, titles, strict=True)):\n",
    "    for fine, exact, values, flux in panels:\n",
    "        value = (exact, values, values - exact, flux)[index]\n",
    "        limit = (common, common, error_max, flux_max)[index]\n",
    "        artist = ax.tripcolor(*fine.points.T, fine.cells, facecolors=value,\n",
    "                             vmin=0.0 if index == 3 else -limit,\n",
    "                             vmax=limit, cmap=\"viridis\" if index == 3 else \"coolwarm\")\n",
    "    for edge in macro.faces:\n",
    "        ax.plot(*macro.points[edge].T, color=\"black\", lw=0.5, alpha=0.7)\n",
    "    ax.set(title=title, xlabel=\"x\", ylabel=\"y\", aspect=\"equal\")\n",
    "    fig.colorbar(artist, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ad6236",
   "metadata": {},
   "source": [
    "## Independently refined classical primal Galerkin baseline\n",
    "\n",
    "The classical baseline uses one globally conforming P2 space on independent 16×16, 32×32 and 64×64 triangular grids, with the same physical coefficient, source and boundary values:\n",
    "\n",
    "$$\n",
    "(\\nabla p_h^{\\mathrm{CG}},\\nabla v_h)=(f,v_h),\\qquad\n",
    "p_h^{\\mathrm{CG}}\\big\\vert_{\\partial\\Omega}=p_D.\n",
    "$$\n",
    "\n",
    "The shared volume operator integrates this already stated form. The explicit essential lifting $A_{ff}p_f=f_f-A_{fb}p_D$ selects its classical boundary values. Check this reference's own errors before comparing fields. The multiscale solve above retains its simpler macro mesh and separate local spaces; it does not use this refined global grid. An exactly representable polynomial has roundoff-level classical errors on all three grids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efbb139c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from threadpoolctl import threadpool_limits\n",
    "from pymhm.fem.scalar.triangle import scalar_operators, nodal_space, tabulate\n",
    "from pymhm.fem.scalar.operators import triangle_quadrature\n",
    "from pymhm.linalg.linear import solve_linear\n",
    "\n",
    "classical_rows = []\n",
    "bary, weights = triangle_quadrature(12)\n",
    "for resolution in (16, 32, 64):\n",
    "    reference_mesh = TriangleMesh.unit_square(resolution)\n",
    "    A, _, forcing = scalar_operators(reference_mesh, 2, diffusion=1.0, source=source, order=12)\n",
    "    dofs, nodes = nodal_space(reference_mesh, 2)\n",
    "    boundary = np.flatnonzero(np.any(np.isclose(nodes, 0.0) | np.isclose(nodes, 1.0), axis=1))\n",
    "    free = np.setdiff1d(np.arange(len(nodes)), boundary)\n",
    "    reference_pressure = np.zeros(len(nodes))\n",
    "    reference_pressure[boundary] = exact_pressure(nodes[boundary])\n",
    "    rhs = forcing[free] - A[free][:, boundary] @ reference_pressure[boundary]\n",
    "    with threadpool_limits(1):\n",
    "        reference_pressure[free] = solve_linear(A[free][:, free], rhs)\n",
    "    _, _, basis, gradient, _ = tabulate(reference_mesh, 2, bary)\n",
    "    points = np.einsum(\"qi,tia->tqa\", bary, reference_mesh.points[reference_mesh.cells])\n",
    "    difference = reference_pressure[dofs] @ basis.T - exact_pressure(points.reshape(-1, 2)).reshape(points.shape[:2])\n",
    "    error = float(np.sqrt(reference_mesh.areas @ (difference**2 @ weights)))\n",
    "    classical_rows.append({\"resolution\": resolution, \"elements\": len(reference_mesh.cells),\n",
    "                           \"pressure_L2\": error})\n",
    "for row in classical_rows:\n",
    "    print(row)\n",
    "# Polynomial patches can be exactly represented, in which case all levels are\n",
    "# at roundoff. Otherwise verify the classical reference's own refinement.\n",
    "assert (max(row[\"pressure_L2\"] for row in classical_rows) < 1e-10\n",
    "        or all(b[\"pressure_L2\"] < a[\"pressure_L2\"] for a, b in zip(classical_rows, classical_rows[1:])))\n",
    "centers = reference_mesh.points[reference_mesh.cells].mean(axis=1)\n",
    "_, _, center_basis, center_gradient, _ = tabulate(reference_mesh, 2, np.full((1, 3), 1/3))\n",
    "reference_values = (reference_pressure[dofs] @ center_basis.T)[:, 0]\n",
    "reference_flux = -np.einsum(\"ti,tqia->tqa\", reference_pressure[dofs], center_gradient)[:, 0]\n",
    "fig, axes = plt.subplots(1, 2, figsize=(9, 4), layout=\"constrained\")\n",
    "for ax, values, title in zip(axes, (reference_values, np.linalg.norm(reference_flux, axis=1)),\n",
    "                             (\"Fine classical P2 pressure\", \"Fine classical Darcy flux magnitude\"), strict=True):\n",
    "    artist = ax.tripcolor(*reference_mesh.points.T, reference_mesh.cells, facecolors=values, cmap=\"viridis\")\n",
    "    for edge in macro.faces:\n",
    "        ax.plot(*macro.points[edge].T, color=\"black\", lw=0.5, alpha=0.7)\n",
    "    ax.set(title=title, xlabel=\"x\", ylabel=\"y\", aspect=\"equal\")\n",
    "    fig.colorbar(artist, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d40ab6",
   "metadata": {},
   "source": [
    "## 7. Refine the independent trace scale\n",
    "\n",
    "Keep the **same 16 macrotriangles** and refine the independent P1 trace segmentation through 1, 2, 4, 8, 16 and 32 pieces. The trace scale is $h_\\Gamma=0.5/s$ for $s$ segments; it is not the fixed macro diameter. The attributed qualification uses P8 local Galerkin pressure with 32 subdivisions and a smooth analytical solution. Its record identifies the executed local fields and the independently checked local-discretization and error-integration controls.\n",
    "\n",
    "The smooth trace comparison in section 6.1 of [Chaumont-Frelet, Paredes and Valentin (2026)](https://doi.org/10.1016/j.camwa.2026.01.016) uses the $\\ell+3/2$ trace estimate. With $\\ell=1$, the target is 2.5 when the finite local Galerkin error is sufficiently smaller than the trace contribution. This is a fixed-macro trace study; the target is not transferred to macro refinement, singular physical domains or an unresolved local solver."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ca2ecf",
   "metadata": {},
   "source": [
    "### Recompute every measured order\n",
    "\n",
    "For successive physical errors $E_{i-1},E_i$ at the actual refinement sizes, compute\n",
    "\n",
    "$$\n",
    "r_i=\\frac{\\log(E_{i-1}/E_i)}{\\log(h_{i-1}/h_i)}.\n",
    "$$\n",
    "\n",
    "The first code lines read one attributed numerical record. They preserve every level and display the complete error/order table. The plotting helper only measures and plots these errors: it constructs no local or global PDE. Targets below are declared from the stated hypotheses, rather than fitted from the data. The final four measured levels give three consecutive orders and a fitted terminal slope. For each declared target $q$, a nearly constant $E/h^q$ provides a second view of the asymptotic regime.\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "820291c0",
   "metadata": {},
   "source": [
    "import json\n",
    "from IPython.display import Markdown, display\n",
    "from examples.tutorial_convergence import refinement_series, asymptotic_summary, plot_method_series\n",
    "refinement_record = ROOT / 'examples/results/unfitted/convergence/rate-verification.json'\n",
    "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n",
    "refinement_family = next((refinement_row for refinement_row in refinement_data['families'] if refinement_row['trace_degree'] == 1))\n",
    "refinement_rows = [dict(H=0.5 / refinement_segments, error=refinement_error) for refinement_segments, refinement_error in zip(refinement_family['segments'], refinement_family['errors'], strict=True)]\n",
    "refinement_error_keys = {'broken gradient L2': 'error'}\n",
    "refinement_targets = {'broken gradient L2': 2.5}\n",
    "refinement_rows = sorted(refinement_rows, key=lambda refinement_row: refinement_row['H'], reverse=True)\n",
    "refinement_sizes = np.asarray([refinement_row['H'] for refinement_row in refinement_rows])\n",
    "refinement_field_errors = {refinement_field: np.asarray([refinement_row[refinement_key] for refinement_row in refinement_rows]) for refinement_field, refinement_key in refinement_error_keys.items()}\n",
    "refinement_orders = {refinement_field: np.log(refinement_values[:-1] / refinement_values[1:]) / np.log(refinement_sizes[:-1] / refinement_sizes[1:]) for refinement_field, refinement_values in refinement_field_errors.items()}\n",
    "refinement_header = ['Refinement size'] + [refinement_column for refinement_field in refinement_error_keys for refinement_column in (refinement_field, 'Order')]\n",
    "refinement_lines = [' | '.join(refinement_header), ' | '.join(['---'] * len(refinement_header))]\n",
    "for refinement_index, refinement_size in enumerate(refinement_sizes):\n",
    "    refinement_values = [f'{refinement_size:.6g}']\n",
    "    for refinement_field in refinement_error_keys:\n",
    "        refinement_values.extend([f'{refinement_field_errors[refinement_field][refinement_index]:.6e}', '—' if refinement_index == 0 else f'{refinement_orders[refinement_field][refinement_index - 1]:.3f}'])\n",
    "    refinement_lines.append(' | '.join(refinement_values))\n",
    "display(Markdown('\\n'.join(refinement_lines)))\n",
    "refinement_series_data = refinement_series(refinement_record, refinement_rows, 'H', refinement_error_keys, refinement_targets, method='unfitted', spaces='P8/r32 locals; segmented P1 traces; 16 fixed macrotriangles', rate_provenance='Chaumont-Frelet, Paredes and Valentin (2026), Theorem 2, equation (20): fixed-macro trace contribution of order 5/2 for q=ell=1 with exact local solves, pressure piecewise H^4 and A piecewise W^{2,infinity} on Lipschitz material regions, A grad(p) in H(div), fitted skeletal subfaces and regular quasi-uniform macrocells; finite P8/r32 local error is controlled separately.', refinement='trace segment size on a fixed macro mesh', root=ROOT)\n",
    "print({'record': refinement_series_data['record'], 'sha256': refinement_series_data['sha256']})\n",
    "for refinement_field, refinement_summary in asymptotic_summary(refinement_series_data).items():\n",
    "    print(refinement_field, {'last_three_orders': [round(refinement_order, 3) for refinement_order in refinement_summary['orders']], 'fitted_order': round(refinement_summary['fitted_order'], 3), 'target': refinement_summary['target'], 'normalized_amplitude_ratio': refinement_summary.get('amplitude_ratio')})\n",
    "plot_method_series(refinement_series_data)\n",
    "plt.show()\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "8ecdeff5",
   "metadata": {},
   "source": [
    "The upper panel preserves all coarse and fine measurements. Dashed lines show declared target powers anchored at the finest measured error. The shaded interval always contains the final four levels; it does not select points to improve a fitted slope. Read the consecutive orders together with the target-normalized errors and the stated quadrature/local-equation controls. A slope alone does not establish the hypotheses of an error estimate."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}