{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "af4a9e2001e1",
   "source": [
    "# Complex Helmholtz MHM: a resolved plane wave and variational workflow\n",
    "\n",
    "Define outgoing physical data, local Helmholtz forms, a normal-flux interface\n",
    "and global pressure continuity by hand. Use generic mesh, form, assembly and\n",
    "solve abstractions; the refinement plots refer to this same declared family.\n",
    "The locked `introduction` environment executes the native UFL calculations.\n",
    "\n",
    "## 1. Derive the wave and its boundary data\n",
    "\n",
    "On the unit square with unit density and modulus, choose $\\omega=10\\pi$,\n",
    "$d=(\\cos(\\pi/13),\\sin(\\pi/13))$ and $p_*(x)=\\exp(i\\omega d\\cdot x)$.\n",
    "The volume source is zero because $\\lVert d\\rVert=1$:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "-\\Delta p-\\omega^2p&=0,\\\\\n",
    "\\partial_np-i\\omega p&=g,\n",
    "&g&=i\\omega(d\\cdot n-1)p_* \\quad\\text{on }\\partial\\Omega.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Represent complex pressure by two interleaved real components $z=(\\Re p,\\Im p)$.\n",
    "The two components are algebraic coordinates, rather than a physical vector.\n",
    "The weak impedance contribution is $\\omega(z_1v_0-z_0v_1)$. This exact real/imaginary\n",
    "embedding keeps the outgoing sign. Exterior data act only on the physical\n",
    "boundary; local interior faces still carry the unknown negative normal gradient.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "519b2848",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import os\n",
    "import sys\n",
    "from pymhm.io.workspace import workspace_from_archive\n",
    "\n",
    "# Download verified support files; this operation does not execute them.\n",
    "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/4ac46f0022294400c5ffc8246fa3bcf42924e7f8dedd6ced12f57bfece4b0bcf/introductory_methods-companion.zip\"\n",
    "COMPANION_SHA256 = \"4ac46f0022294400c5ffc8246fa3bcf42924e7f8dedd6ced12f57bfece4b0bcf\"\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",
    "# Explicitly prepare declared data with the downloaded Python helpers.\n",
    "from scripts.notebook_reproduction import notebook_workspace\n",
    "\n",
    "ROOT = notebook_workspace(\"waves/helmholtz/introductory_methods.ipynb\", directory=ROOT)\n",
    "root = ROOT\n",
    "print(\"Workspace:\", ROOT)\n",
    "\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "81bfd1ea693f",
   "source": [
    "## 2. Write the UFL local forms and trace pairings\n",
    "\n",
    "Choose continuous local Q4 pressure on a $2\\times2$ subdivision of every\n",
    "macro square, and vector P2 coordinates for the complex normal-flux trace.\n",
    "The definition cell below writes the common volume and impedance form, source,\n",
    "local trial pairing and independent global test pairing:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " a_K(p,v)&=(\\nabla p,\\nabla v)_K-\\omega^2(p,v)_K\n",
    "                   -\\langle i\\omega p,v\\rangle_{\\partial K\\cap\\partial\\Omega},\\\\\n",
    " b_K(\\lambda,v)&=\\langle\\lambda,v\\rangle_{\\partial K\\setminus\\partial\\Omega},\\\\\n",
    " c_K(p,\\mu)&=-\\langle p,\\mu\\rangle_{\\partial K\\setminus\\partial\\Omega}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "`local.native_space` owns native coefficient ordering. `local.trace_pairings`\n",
    "owns geometric incidence signs. Declare the mathematical signs directly;\n",
    "do not implement face orientations in the notebook. There are no retained\n",
    "static modes at the selected nonresonant frequency. Positive frequency alone\n",
    "does not establish local invertibility at arbitrary frequencies.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "id": "32a7fb10003e",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import ufl\n",
    "import basix\n",
    "import basix.ufl\n",
    "from types import SimpleNamespace\n",
    "from pymhm import (Equation, TriangleMesh, FaceSpace, SkeletonSpace, MeshHierarchy,\n",
    "                   LocalContext, bind_interface, bind_problem, assemble, solve,\n",
    "                   MultiscaleProblem, compile_form)\n",
    "from pymhm.meshes.cartesian import CartesianMacroMesh\n",
    "from pymhm.backends.spaces import bind_space, coefficient_map\n",
    "from pymhm.fem.scalar.quadrilateral import qk_space, quadrilateral_quadrature\n",
    "from examples.minimal_flow_originals import original_diagnostics\n",
    "\n",
    "wave_frequency=10*np.pi\n",
    "wave_direction=np.array((np.cos(np.pi/13),np.sin(np.pi/13)))\n",
    "\n",
    "def wave_truth(points):\n",
    "    \"\"\"Evaluate the independently specified unit-amplitude plane wave as Re/Im.\"\"\"\n",
    "    phase=wave_frequency*(np.asarray(points)@wave_direction)\n",
    "    return np.column_stack((np.cos(phase),np.sin(phase)))\n",
    "\n",
    "def wave_gradient(points):\n",
    "    \"\"\"Differentiate the analytical wave in physical Cartesian coordinates.\"\"\"\n",
    "    phase=wave_frequency*(np.asarray(points)@wave_direction)\n",
    "    return wave_frequency*np.stack((-np.sin(phase),np.cos(phase)),axis=1)[:,:,None]*wave_direction\n",
    "\n",
    "def wave_forms(V,domain):\n",
    "    \"\"\"Write the Helmholtz volume and outgoing physical-boundary impedance forms.\"\"\"\n",
    "    z,v=ufl.TrialFunction(V),ufl.TestFunction(V)\n",
    "    x=ufl.SpatialCoordinate(domain)\n",
    "    n=ufl.FacetNormal(domain)\n",
    "    exterior=ufl.conditional(ufl.Or(ufl.Or(ufl.lt(x[0],1e-12),ufl.gt(x[0],1-1e-12)),\n",
    "                                  ufl.Or(ufl.lt(x[1],1e-12),ufl.gt(x[1],1-1e-12))),1.,0.)\n",
    "    phase=wave_frequency*(x[0]*wave_direction[0]+x[1]*wave_direction[1])\n",
    "    normal_factor=wave_frequency*(ufl.dot(ufl.as_vector(wave_direction),n)-1)\n",
    "    boundary_source=ufl.as_vector((-normal_factor*ufl.sin(phase),normal_factor*ufl.cos(phase)))\n",
    "    dx=ufl.Measure('dx',domain=domain,metadata={'quadrature_degree':16})\n",
    "    ds=ufl.Measure('ds',domain=domain,metadata={'quadrature_degree':16})\n",
    "    a=(ufl.inner(ufl.grad(z),ufl.grad(v))-wave_frequency**2*ufl.inner(z,v))*dx\n",
    "    a+=wave_frequency*(z[1]*v[0]-z[0]*v[1])*exterior*ds\n",
    "    L=ufl.inner(boundary_source,v)*exterior*ds\n",
    "    return a,L,exterior\n",
    "\n",
    "wave_element=basix.ufl.element('Lagrange','quadrilateral',4,shape=(2,),\n",
    "                               lagrange_variant=basix.LagrangeVariant.equispaced)\n",
    "\n",
    "def wave_local(local: LocalContext):\n",
    "    \"\"\"Declare Q4 complex pressure, P2 normal moments and outgoing exterior data.\"\"\"\n",
    "    native=local.native_space(wave_element)\n",
    "    z,v=ufl.TrialFunction(native.space),ufl.TestFunction(native.space)\n",
    "    a,L,exterior=wave_forms(native.space,native.mesh)\n",
    "    b=local.trace_pairings(lambda phi,ds:ufl.inner(phi,v)*(1-exterior)*ds)\n",
    "    c=local.trace_pairings(lambda phi,ds:-ufl.inner(phi,z)*(1-exterior)*ds,axis='rows')\n",
    "    local.field('pressure',native)\n",
    "    return local.equations(a=a,L=L,b=b,c=c)\n",
    "\n",
    "def wave_errors(fields,order):\n",
    "    \"\"\"Integrate physical relative pressure/gradient errors with independent quadrature.\"\"\"\n",
    "    reference,weights=quadrilateral_quadrature(order)\n",
    "    squared=np.zeros(2)\n",
    "    for field in fields:\n",
    "        fine=field.mesh\n",
    "        origins=fine.points[fine.cells[:,0]]\n",
    "        points=(origins[:,None,:]+reference[None,:,:]*fine.spacing).reshape(-1,2)\n",
    "        cell_ids=np.repeat(np.arange(len(fine.cells)),len(reference))\n",
    "        values,gradient=field.values_and_gradient(points,cells=cell_ids)\n",
    "        physical_weights=(fine.areas[:,None]*weights[None,:]).ravel()\n",
    "        squared[0]+=np.sum(physical_weights*np.sum((values-wave_truth(points))**2,axis=1))\n",
    "        squared[1]+=np.sum(physical_weights*np.sum((gradient-wave_gradient(points))**2,axis=(1,2)))\n",
    "    return {'pressure_relative_l2':float(np.sqrt(squared[0])),\n",
    "            'gradient_relative_l2':float(np.sqrt(squared[1])/wave_frequency)}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "315ecf3f438d",
   "source": [
    "## 3. Define the meshes and global equation, assemble and solve\n",
    "\n",
    "The multiplier is the negative physical normal gradient on interior faces.\n",
    "Exterior multiplier coordinates are zero placeholders because the Robin data\n",
    "are already in the local impedance form and load. They are not an assertion\n",
    "that the physical exterior flux vanishes. Fix only these placeholders to zero.\n",
    "The global equation sums incident pressure moments with the negative sign\n",
    "chosen in $c_K$; its right-hand side is zero.\n",
    "\n",
    "After `assemble` and `solve`, evaluate the named pressure fields using their\n",
    "executed bases. Independent quadrature orders 10 and 14 integrate physical\n",
    "relative pressure and broken gradient errors. The separate original-row\n",
    "checker measures uncondensed coefficient equations, not a physical field norm.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "id": "1f150a2748d2",
   "execution_count": null,
   "outputs": [],
   "source": [
    "wave_macro=CartesianMacroMesh(12)\n",
    "wave_skeleton=SkeletonSpace(wave_macro,tuple(FaceSpace.uniform(2) for _ in wave_macro.faces),components=2)\n",
    "wave_hierarchy=MeshHierarchy(wave_macro,tuple(wave_macro.submesh(i,2) for i in range(len(wave_macro.cells))))\n",
    "wave_interface=bind_interface(wave_skeleton,convention='normal')\n",
    "wave_problem=bind_problem(wave_hierarchy,wave_interface,wave_local,global_equation=Equation(0,0))\n",
    "wave_system=assemble(wave_problem)\n",
    "wave_fixed={int(dof):0. for face in wave_macro.boundary_faces for dof in wave_skeleton.dofs(int(face))}\n",
    "wave_solution=solve(wave_system,fixed=wave_fixed)\n",
    "wave_pressure=wave_solution.field('pressure')\n",
    "wave_norms=wave_errors(wave_pressure,10)\n",
    "wave_norms_finer=wave_errors(wave_pressure,14)\n",
    "wave_original=original_diagnostics(SimpleNamespace(hybrid=wave_solution),\n",
    "    SimpleNamespace(system=wave_system,applied_boundary=-wave_system.global_load[:wave_system.trace_size],\n",
    "                    fixed=wave_fixed,mean_weights=[],mean_values=[],constraints=[]),\n",
    "    [{'complex_pressure_native_block':len(values)} for values in wave_solution.fields])\n",
    "assert wave_original['accepted']\n",
    "assert all(np.isclose(wave_norms[k],wave_norms_finer[k],rtol=1e-7,atol=0) for k in wave_norms)\n",
    "print('By-hand UFL plane-wave errors:',wave_norms,flush=True)\n",
    "print('Original local/global rows:',{k:v for k,v in wave_original.items() if k!='blocks'},flush=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "c3991f9add46",
   "source": [
    "## 4. Independently refine a classical Galerkin reference\n",
    "\n",
    "Use the same `wave_forms` on one globally conforming Q4 mesh, with the same\n",
    "outgoing physical data and no MHM condensation. `MultiscaleProblem.from_global`\n",
    "allows the same general assembly/solve workflow. The native adapter closes its\n",
    "resources explicitly; the portable coefficient map and declared field preserve\n",
    "the executed nodal basis for evaluation.\n",
    "\n",
    "Refine the reference on global grids 16, 32 and 64 before plotting its finest\n",
    "field. Both analytical norms must decrease; a fine reference remains numerical.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "id": "254b6cdac901",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def classical_wave(n):\n",
    "    \"\"\"Independently assemble conforming Q4 Galerkin with the same outgoing UFL data.\"\"\"\n",
    "    fine=CartesianMacroMesh(n)\n",
    "    native=bind_space(fine,element=wave_element)\n",
    "    try:\n",
    "        a,L,_=wave_forms(native.space,native.mesh)\n",
    "        problem=MultiscaleProblem.from_global(Equation(a,L),native.size)\n",
    "        system=assemble(problem)\n",
    "        state=solve(system)\n",
    "        _,points=qk_space(fine,4)\n",
    "        canonical=state.trace[coefficient_map(native)]\n",
    "    finally:\n",
    "        native.close()\n",
    "    from pymhm.postprocessing.nodal import nodal_field\n",
    "    from pymhm.postprocessing.fields import DiscreteField\n",
    "    field=DiscreteField(nodal_field('pressure',fine,4,components=2),canonical)\n",
    "    errors=wave_errors((field,),10)\n",
    "    errors_finer=wave_errors((field,),14)\n",
    "    assert all(np.isclose(errors[k],errors_finer[k],rtol=1e-7,atol=0) for k in errors)\n",
    "    assert state.raw_residual<=1e-10\n",
    "    return field,{'n':n,**errors,'global_compatibility_relative':state.raw_residual}\n",
    "\n",
    "wave_references=[]\n",
    "for n in (16,32,64):\n",
    "    wave_reference,row=classical_wave(n)\n",
    "    wave_references.append(row)\n",
    "assert all(wave_references[i+1][name]<wave_references[i][name]\n",
    "           for i in range(2) for name in wave_norms)\n",
    "print('Independent global Q4 references:',wave_references,flush=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "169163fd701a",
   "source": [
    "## 5. Plot the complex wave and its errors\n",
    "\n",
    "Each macrocell is sampled independently. The actual macro mesh is drawn on\n",
    "analytical, multiscale, reference and error panels. Individual colorbars remain\n",
    "independent; reconstructed interface values are not averaged.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "id": "403526943256",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from examples.introduction.vector import plot_field_panels\n",
    "from pathlib import Path\n",
    "display=CartesianMacroMesh(8)\n",
    "display_triangles=np.vstack((display.cells[:,(0,1,2)],display.cells[:,(0,2,3)]))\n",
    "points=[];triangles=[];pressure=[];offset=0\n",
    "for cell,field in enumerate(wave_pressure):\n",
    "    vertices=wave_macro.points[wave_macro.cells[cell]]\n",
    "    mapped=vertices[0]+display.points*(vertices[2]-vertices[0])\n",
    "    points.append(mapped);triangles.append(display_triangles+offset)\n",
    "    pressure.append(field.evaluate(mapped));offset+=len(mapped)\n",
    "points=np.concatenate(points);triangles=np.concatenate(triangles);pressure=np.concatenate(pressure)\n",
    "truth=wave_truth(points);reference=wave_reference.evaluate(points)\n",
    "fig=plot_field_panels(wave_macro,{\n",
    "    'Analytical Re pressure':(points,triangles,truth[:,0]),\n",
    "    'MHM Re pressure':(points,triangles,pressure[:,0]),\n",
    "    'Fine Galerkin Re pressure':(points,triangles,reference[:,0]),\n",
    "    'Analytical Im pressure':(points,triangles,truth[:,1]),\n",
    "    'MHM pressure error magnitude':(points,triangles,np.linalg.norm(pressure-truth,axis=1)),\n",
    "    'Galerkin pressure error magnitude':(points,triangles,np.linalg.norm(reference-truth,axis=1)),\n",
    "},figsize=(12,8))\n",
    "output=Path('build/tutorial-methods');output.mkdir(parents=True,exist_ok=True)\n",
    "fig.savefig(output/'helmholtz-wave-fields.png',dpi=200,\n",
    "    metadata={'Author':'IPES Research Group','License':'CC BY 4.0 — https://creativecommons.org/licenses/by/4.0/'})\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "472c4f0d1a80",
   "source": [
    "## 6. Reach the fixed-frequency asymptotic regime\n",
    "\n",
    "The recorded family uses exactly the plane wave, Robin convention, square macro\n",
    "meshes, local Q4/refinement two and polynomial P2 traces defined above. It\n",
    "retains all $n=12,16,24,32,48,64$ observations. The native n12 solve above also\n",
    "checks that its two physical norms agree with the attributed record.\n",
    "\n",
    "For the smooth fixed-frequency regime, Chaumont-Frelet and Valentin (2020),\n",
    "Section 6.3, Figures 6.6 and 6.7, report pressure order $\\ell+2$ and broken\n",
    "gradient order $\\ell+1$ for polynomial traces. Here $\\ell=2$: target orders are\n",
    "four and three. The final three pressure orders are 3.9910, 3.9918 and 3.9930;\n",
    "gradient orders are 2.9906, 2.9928 and 2.9941. The final four-level normalized\n",
    "amplitudes vary by less than one percent.\n",
    "\n",
    "The record distinguishes PyMHM analytical errors from values digitized from\n",
    "the published figures. Local degrees and refinements are explicitly declared\n",
    "here; equality with an unidentified historical local mesh is not claimed.\n",
    "These rates are not a frequency-uniform pollution or heterogeneous resonance\n",
    "estimate.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "id": "2f47ea9deba8",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import json\n",
    "from examples.tutorial_convergence import refinement_series,plot_method_series,asymptotic_summary\n",
    "record_path=ROOT/'examples/results/helmholtz-article/published-convergence.json'\n",
    "record=json.loads(record_path.read_text())\n",
    "selected=[row for row in record['rows'] if row['ell']==2 and row['basis']=='polynomial']\n",
    "rows=[]\n",
    "for resolution in sorted({row['n'] for row in selected}):\n",
    "    fields={row['field']:row['pymhm_relative_error'] for row in selected if row['n']==resolution}\n",
    "    rows.append({'n':resolution,'H':1/resolution,**fields})\n",
    "first=next(row for row in rows if row['n']==12)\n",
    "for name,field in [('pressure_relative_l2','pressure'),('gradient_relative_l2','gradient')]:\n",
    "    np.testing.assert_allclose(wave_norms[name],first[field],rtol=1e-8,atol=0)\n",
    "series=refinement_series(record_path,rows,'H',\n",
    "    {'relative pressure L2':'pressure','relative broken gradient L2':'gradient'},\n",
    "    {'relative pressure L2':4.0,'relative broken gradient L2':3.0},method='helmholtz',\n",
    "    spaces='Q4/r2 local pressure; polynomial P2 complex normal traces; all exterior outgoing Robin data',\n",
    "    rate_provenance='Chaumont-Frelet and Valentin (2020), theorem 5.4: conditional fixed-frequency polynomial one-level rates require global uniqueness, aligned smooth material and frequency-dependent mesh smallness; finite local errors are separate; PyMHM measurements, not digitized article values',\n",
    "    refinement='macro grid spacing',root=ROOT)\n",
    "print('Attributed input:',series['record'],series['sha256'])\n",
    "print(asymptotic_summary(series))\n",
    "figure=plot_method_series(series)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "e5cf33f6d109",
   "source": [
    "## Further control: an exactly represented affine field\n",
    "\n",
    "The following low-frequency complex P2 patch preserves the earlier polynomial\n",
    "reproduction example. Its geometric provider is a convenience adapter for\n",
    "explicit volume and boundary forms. It has no meaningful fitted convergence\n",
    "rate and is distinct from the resolved plane-wave experiment above.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a817cb2",
   "metadata": {},
   "source": [
    "## Declare and solve complex Helmholtz equations\n",
    "\n",
    "The primary native forms use the exact real/imaginary embedding of complex\n",
    "pressure; no complex coefficient is discarded. Unit density/modulus give\n",
    "$-\\Delta p-\\omega^2p=f$ with $\\omega=1/2$ and the independent affine\n",
    "source $f=-\\omega^2p$. This chosen finite patch has invertible local operators;\n",
    "positive frequency alone does not exclude local resonances in other inputs.\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_T(p,v)&=(\\nabla p,\\nabla v)_T-\\omega^2(p,v)_T,\\\\\n",
    "B_T\\lambda&=\\langle\\lambda,v\\rangle_{\\partial T},\\\\\n",
    "C_Tp&=-\\langle p,\\mu\\rangle_{\\partial T}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The signed trace is $-\\nabla p\\cdot n$ in two interleaved real components.\n",
    "The global RHS explicitly prescribes negative complex Dirichlet moments.\n",
    "The native example has no absorption, PML or resonance-avoidance claim.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9552127b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib.util\n",
    "from pymhm import Equation, LocalEquations, MultiscaleProblem, assemble, columns, rows, compile_form\n",
    "from pymhm.meshes.triangle import TriangleMesh\n",
    "from pymhm.fem.traces.interval import FaceSpace, SkeletonSpace\n",
    "from pymhm.fem.scalar.operators import boundary_data\n",
    "from examples.variational_darcy import NativeFormProvider\n",
    "\n",
    "native_available = importlib.util.find_spec(\"dolfinx\") is not None\n",
    "if native_available:\n",
    "    import ufl\n",
    "    from dolfinx import fem\n",
    "\n",
    "    mesh = TriangleMesh.unit_square()\n",
    "    trace_space = SkeletonSpace(mesh, components=2)\n",
    "    omega = 0.5\n",
    "\n",
    "    def exact_components(points):\n",
    "        \"\"\"Return real/imaginary components of (1+i)(1+x+2y).\"\"\"\n",
    "        scalar = 1 + points[:, 0] + 2 * points[:, 1]\n",
    "        return np.column_stack((scalar, scalar))\n",
    "\n",
    "    def local_forms(context):\n",
    "        \"\"\"Write both real components of the complex Helmholtz equation.\"\"\"\n",
    "        space = fem.functionspace(context.space.mesh, (\"Lagrange\", 2, (2,)))\n",
    "        u, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n",
    "        x = ufl.SpatialCoordinate(context.space.mesh)\n",
    "        dx = ufl.dx(domain=context.space.mesh)\n",
    "        exact = ufl.as_vector((1 + x[0] + 2 * x[1], 1 + x[0] + 2 * x[1]))\n",
    "        return LocalEquations(\n",
    "            (ufl.inner(ufl.grad(u), ufl.grad(v)) - omega**2 * ufl.inner(u, v)) * dx,\n",
    "            ufl.inner(-(omega**2) * exact, v) * dx,\n",
    "            columns(\n",
    "                *(\n",
    "                    float(sign) * v[j] * context.ds(side + 1)\n",
    "                    for side, sign in enumerate(context.signs)\n",
    "                    for j in range(2)\n",
    "                )\n",
    "            ),\n",
    "            rows(\n",
    "                *(\n",
    "                    -float(sign) * u[j] * context.ds(side + 1)\n",
    "                    for side, sign in enumerate(context.signs)\n",
    "                    for j in range(2)\n",
    "                )\n",
    "            ),\n",
    "            context.trace_dofs,\n",
    "            metadata={\"points\": space.tabulate_dof_coordinates()[:, :2].copy()},\n",
    "        )\n",
    "\n",
    "    boundary, _ = boundary_data(trace_space, exact_components)\n",
    "    problem = MultiscaleProblem(\n",
    "        Equation(0, -boundary),\n",
    "        NativeFormProvider(mesh, trace_space, local_forms, 2),\n",
    "        range(len(mesh.cells)),\n",
    "        trace_space.size,\n",
    "        (0,) * len(mesh.cells),\n",
    "    )\n",
    "    system = assemble(problem)\n",
    "    solution = system.solve()\n",
    "    errors = [\n",
    "        float(np.max(np.abs(field.reshape(-1, 2) - exact_components(record[\"points\"]))))\n",
    "        for field, record in zip(solution.fields, system.local_metadata, strict=True)\n",
    "    ]\n",
    "    assert max(errors) < 1e-10\n",
    "    print(\n",
    "        {\n",
    "            \"native_available\": True,\n",
    "            \"spaces\": \"complex P2 pressure / P0 complex trace\",\n",
    "            \"maximum_component_nodal_error_by_cell\": errors,\n",
    "            \"global_original_relative_residual\": solution.raw_residual,\n",
    "        }\n",
    "    )\n",
    "else:\n",
    "    print(\n",
    "        {\n",
    "            \"native_available\": False,\n",
    "            \"status\": \"User-written UFL primary solve not executed: select the Pixi fem kernel.\",\n",
    "        }\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "649dcc3f",
   "metadata": {},
   "source": [
    "## Compare the supported methods\n",
    "These comparisons invoke the existing formulation implementations through the application catalogue and preserve their numerical records. The user-defined primary UFL solve above exercises the generic local/global API; the catalogue drivers remain compatibility comparisons.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46163fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from examples.tutorial_scalar_variants import (\n",
    "    VARIANTS,\n",
    "    affine_pressure,\n",
    "    affine_flux,\n",
    "    run_variant,\n",
    ")\n",
    "\n",
    "record = run_variant(\"helmholtz\")\n",
    "print(record[\"spaces_and_conventions\"])\n",
    "print(record[\"metrics\"])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "physical-plot-description",
   "metadata": {},
   "source": [
    "## Plot pressure components and error with the macro mesh\n",
    "\n",
    "Both components of this affine complex pressure are identical. Sampled nodal values reproduce the affine field; the panel triangulation displays those values separately in each macrocell. The macro mesh is drawn on every panel.\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "physical-plot-code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "assert native_available, \"Select the locked introduction environment for this UFL demonstration.\"\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.tri import Triangulation\n",
    "\n",
    "fig, axes = plt.subplots(2, 3, figsize=(10, 6), layout=\"constrained\")\n",
    "for component in range(2):\n",
    "    for column in range(3):\n",
    "        ax = axes[component, column]\n",
    "        for values, record in zip(solution.fields, system.local_metadata, strict=True):\n",
    "            points = record[\"points\"]\n",
    "            exact = exact_components(points)[:, component]\n",
    "            numerical = values.reshape(-1, 2)[:, component]\n",
    "            shown = (exact, numerical, np.abs(numerical - exact))[column]\n",
    "            image = ax.tripcolor(Triangulation(points[:, 0], points[:, 1]), shown,\n",
    "                                 shading=\"gouraud\", vmin=0 if column==2 else 1,\n",
    "                                 vmax=1e-10 if column==2 else 4)\n",
    "        ax.triplot(mesh.points[:, 0], mesh.points[:, 1], mesh.cells, color=\"black\", linewidth=0.6)\n",
    "        name = \"Real\" if component==0 else \"Imaginary\"\n",
    "        ax.set(title=f\"{name}: {('exact', 'MHM', 'absolute error')[column]}\",\n",
    "               xlabel=\"x\", ylabel=\"y\", aspect=\"equal\")\n",
    "        fig.colorbar(image, ax=ax, label=\"Pressure\" if column<2 else \"Pressure error\")\n",
    "Path(\"build/tutorial-methods\").mkdir(parents=True, exist_ok=True)\n",
    "fig.savefig(\"build/tutorial-methods/helmholtz-affine-fields.png\", dpi=200, metadata={\"Author\":\"IPES Research Group\", \"License\":\"CC BY 4.0 — https://creativecommons.org/licenses/by/4.0/\"})\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "classical-reference-description",
   "metadata": {},
   "source": [
    "## Independent conforming Galerkin reference\n",
    "\n",
    "Assemble the same real/imaginary Helmholtz forms on global P2 grids of 16, 32 and 64 subdivisions. Prescribe the exact exterior pressure strongly and solve with the common sparse factorization owner. The affine field is exactly representable, so refinement verifies reproduction instead of supplying a fitted rate. The MHM example uses only two macrocells.\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "classical-reference-code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "from mpi4py import MPI\n",
    "from dolfinx import mesh as native_mesh\n",
    "from pymhm.linalg.linear import factorize\n",
    "\n",
    "classical = []\n",
    "for resolution in (16, 32, 64):\n",
    "    domain = native_mesh.create_unit_square(MPI.COMM_SELF, resolution, resolution)\n",
    "    space = fem.functionspace(domain, (\"Lagrange\", 2, (2,)))\n",
    "    u, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n",
    "    x = ufl.SpatialCoordinate(domain)\n",
    "    dx = ufl.dx(domain=domain)\n",
    "    exact = ufl.as_vector((1+x[0]+2*x[1], 1+x[0]+2*x[1]))\n",
    "    a = (ufl.inner(ufl.grad(u), ufl.grad(v)) - omega**2*ufl.inner(u,v))*dx\n",
    "    L = -omega**2 * ufl.inner(exact, v)*dx\n",
    "    points = space.tabulate_dof_coordinates()[:, :2]\n",
    "    size = 2*len(points)\n",
    "    matrix, load = compile_form(a, (size,size)), compile_form(L, (size,))\n",
    "    exterior = np.any(np.isclose(points,0) | np.isclose(points,1), axis=1)\n",
    "    fixed = (2*np.flatnonzero(exterior)[:,None] + np.arange(2)).ravel()\n",
    "    free = np.setdiff1d(np.arange(size),fixed)\n",
    "    values = np.zeros(size)\n",
    "    values[fixed] = exact_components(points).ravel()[fixed]\n",
    "    with factorize(matrix[free][:,free]) as factor:\n",
    "        values[free] = factor.solve(load[free] - matrix[free][:,fixed] @ values[fixed])\n",
    "    reference = fem.Function(space)\n",
    "    reference.x.array[:] = values\n",
    "    difference = reference-exact\n",
    "    l2 = np.sqrt(fem.assemble_scalar(fem.form(ufl.inner(difference,difference)*dx)))\n",
    "    assert l2 < 1e-10\n",
    "    classical.append({\"resolution\":resolution,\"pressure_L2\":float(l2)})\n",
    "classical\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "id": "ea4770a5f6fa",
   "source": [
    "## References\n",
    "\n",
    "- Chaumont-Frelet and Valentin (2020), *A Multiscale Hybrid-Mixed Method for the Helmholtz Equation in Heterogeneous Domains*. [DOI:10.1137/19M1255616](https://doi.org/10.1137/19M1255616).\n",
    "\n",
    "The fixed-frequency analytical qualification declares its own local discretization;\n",
    "it does not claim a matched historical table without its original inputs.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "pymhm": {
   "kind": "introductory",
   "methods": [
    "complex hybrid Helmholtz"
   ],
   "problem": "waves/helmholtz",
   "scope": "analytical patches; no convergence or literature-reproduction claim",
   "variants": [
    "helmholtz"
   ]
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}