{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78b4441d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "composition mismatch L2: 1.795339620998114e-05\n",
      "Newton iteration  0\n",
      "err =  2.29107865967357\n",
      "Newton iteration  1\n",
      "err =  0.25871216757154775\n",
      "Newton iteration  2\n",
      "err =  0.042037219808741005\n",
      "Newton iteration  3\n",
      "err =  0.001753044966358795\n",
      "Newton iteration  4\n",
      "err =  3.3228273248850992e-06\n",
      "Newton iteration  5\n",
      "err =  1.1983821811469996e-11\n",
      "Newton iteration  6\n",
      "err =  8.570424917112041e-16\n",
      "Newton iteration  0\n",
      "err =  2.29107865967357\n",
      "Newton iteration  1\n",
      "err =  3.443394500729945e-15\n",
      "symbolic error:  3.794941764432178e-06\n",
      "composed error:  0.2500811637999869\n"
     ]
    }
   ],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.fem import CoordinateTrafo\n",
    "from ngsolve.solvers import Newton\n",
    "from netgen.geom2d import unit_square\n",
    "from netgen.meshing import Mesh as NGMesh, MeshPoint, Element1D, Element0D, Pnt\n",
    "import numpy as np\n",
    "\n",
    "# 1D mesh for the tabulated nonlinearity Pi(s) = s², s in [-2, 2]\n",
    "n = 200\n",
    "ngmesh1d = NGMesh()\n",
    "ngmesh1d.dim = 1\n",
    "pts = [ngmesh1d.Add(MeshPoint(Pnt(-2 + 4*i/n, 0, 0))) for i in range(n + 1)]\n",
    "for i in range(n):\n",
    "    ngmesh1d.Add(Element1D([pts[i], pts[i + 1]], index=1))\n",
    "ngmesh1d.Add(Element0D(pts[0], index=1))\n",
    "ngmesh1d.Add(Element0D(pts[-1], index=2))\n",
    "\n",
    "mesh1d = Mesh(ngmesh1d)\n",
    "fes1d = H1(mesh1d, order=3)\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "fes = H1(mesh, order=3)\n",
    "\n",
    "pi_gf = GridFunction(fes1d)\n",
    "pi_gf.Set(x**2)\n",
    "\n",
    "# Compoition function computes Pi(u) by evaluating the 1D GridFunction pi_gf at the values of u.\n",
    "def compose_gridfunction(pi, u):\n",
    "    u_clipped = IfPos(u - 2, 2, IfPos(-2 - u, -2, u))\n",
    "    pi_of_u = pi(\n",
    "        CoordinateTrafo(CF((u_clipped, 0, 0)), mesh1d.Materials(\".*\"))\n",
    "    )\n",
    "    result = GridFunction(fes)\n",
    "    result.Set(pi_of_u)\n",
    "    return result\n",
    "\n",
    "u_exact = cos(pi*x) * cos(pi*y)\n",
    "\n",
    "# PDE: -Delta(u) + u + Pi(u) = phi, with Pi(u) = u²\n",
    "phi = (2*pi**2 + 1) * u_exact + u_exact**2\n",
    "\n",
    "# Error caused solely by tabulating/interpolating Pi(s) = s² on mesh1d\n",
    "pi_of_u_exact = compose_gridfunction(pi_gf, u_exact)\n",
    "composition_mismatch = pi_of_u_exact - u_exact**2\n",
    "\n",
    "print(\n",
    "    \"composition mismatch L2:\",\n",
    "    sqrt(Integrate(composition_mismatch**2, mesh)),\n",
    ")\n",
    "\n",
    "u, v = fes.TnT()\n",
    "\n",
    "# Control case: symbolic Pi(u) = u². Native Newton works.\n",
    "a_symbolic = BilinearForm(fes)\n",
    "a_symbolic += (InnerProduct(grad(u), grad(v)) + u*v + u**2*v - phi*v) * dx\n",
    "\n",
    "gfu_symbolic = GridFunction(fes)\n",
    "Newton(a_symbolic, gfu_symbolic, fes.FreeDofs(), maxit=20, maxerr=1e-12)\n",
    "\n",
    "# Problem case: Pi is represented by a 1D GridFunction and composed with u.\n",
    "a_composed = BilinearForm(fes)\n",
    "a_composed += (\n",
    "    InnerProduct(grad(u), grad(v))\n",
    "    + u*v\n",
    "    + compose_gridfunction(pi_gf, u)*v\n",
    "    - phi*v\n",
    ") * dx\n",
    "\n",
    "gfu_composed = GridFunction(fes)\n",
    "Newton(a_composed, gfu_composed, fes.FreeDofs(), maxit=20, maxerr=1e-12)\n",
    "\n",
    "print(\"symbolic error: \",\n",
    "      sqrt(Integrate((gfu_symbolic - u_exact)**2, mesh)))\n",
    "print(\"composed error: \",\n",
    "      sqrt(Integrate((gfu_composed - u_exact)**2, mesh)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_ngsolve",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}