{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# FEEC basics\n",
    "\n",
    "This tutorial provides a practical introduction to working with finite element exterior calculus (FEEC) in Struphy. It covers the fundamental concepts and tools you need to set up and manipulate the discrete de Rham complex.\n",
    "\n",
    "**Prerequisite knowledge**: Familiarity with the mathematical concepts of FEEC is helpful. For comprehensive background, see [the numerics section of the Struphy documentation](https://struphy-hub.github.io/struphy/sections/subsections/numerics-geomFE.html) and the references therein.\n",
    "\n",
    "**What you'll learn**:\n",
    "1. How to set up a `Derham` object and understand its key attributes\n",
    "2. How to access 1D spline space information via `SplineSpace1D`\n",
    "3. How to create callable spline functions for evaluation\n",
    "4. How to use differential operators (grad, curl, div)\n",
    "5. How to work with geometric projectors and the commuting property\n",
    "\n",
    "This tutorial is the foundation for more advanced FEEC work, such as boundary condition handling (see the [boundary conditions tutorial](dev_tutorial_feec_bcs.ipynb)).\n",
    "\n",
    "## 1. Setting up a Derham object\n",
    "\n",
    "The `Derham` class creates the discrete de Rham complex—a fundamental data structure in FEEC that contains all the finite element spaces and operators you need.\n",
    "\n",
    "To create a `Derham` object, you need two things:\n",
    "1. A spatial grid (`TensorProductGrid`)\n",
    "2. Configuration options (`DerhamOptions`)\n",
    "\n",
    "Let's start with a simple 1D domain with 16 elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from struphy.feec.psydac_derham import Derham\n",
    "from struphy.topology.grids import TensorProductGrid\n",
    "from struphy.io.options import DerhamOptions\n",
    "\n",
    "# Create a 1D grid with 16 elements\n",
    "grid = TensorProductGrid(num_elements=(16, 1, 1))\n",
    "\n",
    "# Create default Derham object with default options\n",
    "derham_opts = DerhamOptions()\n",
    "derham = Derham(grid, derham_opts)\n",
    "\n",
    "print(\"Derham object created successfully\")\n",
    "print(f\"Grid has {derham.grid.num_elements} elements\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "### Default configuration\n",
    "\n",
    "By default, the `Derham` object uses spline degree $p=1$ in each direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Spline degrees: {derham.degree}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "And periodic boundary conditions (no constraints) in all directions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Boundary conditions: {derham.bcs}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Let's also check the quadrature configuration used for numerical integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Quadrature points per direction: {derham.nquads}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## 2. Understanding Derham attributes and SplineSpace1D\n",
    "\n",
    "The `Derham` object provides several important attributes for accessing the discrete spaces and their properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "### 2.1 Space identifiers (`space_to_form`)\n",
    "\n",
    "FEEC organizes spaces by the differential form degree they support. The mapping from space identifier to form degree is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Space to form mapping:\")\n",
    "for space_id, form_degree in derham.space_to_form.items():\n",
    "    print(f\"  {space_id:4s} -> form degree {form_degree}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### 2.2 Discrete function spaces (`Vh`)\n",
    "\n",
    "The `Vh` dictionary contains the vector spaces (from Psydac's `psydac.linalg`) where finite element coefficients live. Each space corresponds to one of the de Rham spaces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Discrete spaces (Vh):\")\n",
    "for space_id, V in derham.coeff_spaces.items():\n",
    "    if space_id in derham.space_to_form:\n",
    "        print(f\"  {space_id:4s}: {V}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "### 2.3 Polar spaces (`Vh_pol`)\n",
    "\n",
    "In addition to Cartesian spaces, FEEC can handle polar coordinates. These spaces are available via `Vh_pol` (if they have been initialized):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Polar spaces (Vh_pol):\")\n",
    "if derham.polar_coeff_spaces is not None:\n",
    "    for space_id, V in derham.polar_coeff_spaces.items():\n",
    "        if space_id in derham.space_to_form:\n",
    "            print(f\"  {space_id:4s}: {V}\")\n",
    "else:\n",
    "    print(\"  Polar spaces not initialized (None)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "### 2.4 1D spline space information (`SplineSpace1D`)\n",
    "\n",
    "One powerful feature of Struphy is the ability to access information about the 1D spline spaces underlying the 3D tensor product construction. Each de Rham space corresponds to a different combination of 1D spline space types in the three directions.\n",
    "\n",
    "The spline space types are accessed via attributes like `V0splines`, `V1splines`, `V2splines`, and `V3splines`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access the 0-form (H1) spline space information\n",
    "V0_splines = derham.V0splines\n",
    "print(\"H1 space (0-forms):\")\n",
    "print(f\"  Type: {type(V0_splines)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "These `SplineSpace1D` objects contain detailed information about the 1D tensor product construction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Each direction's spline space info is stored in separate components\n",
    "print(\"\\nQuadrature grid points per direction (for H1 space):\")\n",
    "quad_pts = V0_splines.quad_grid_pts\n",
    "for i, pts_tuple in enumerate(quad_pts):\n",
    "    if isinstance(pts_tuple, tuple):\n",
    "        print(f\"  Direction {i}: {len(pts_tuple)} components\")\n",
    "        for j, pts in enumerate(pts_tuple):\n",
    "            print(f\"    Component {j}: {len(pts)} quadrature points\")\n",
    "    else:\n",
    "        print(f\"  Direction {i}: {len(pts_tuple)} quadrature points\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "Let's examine spline space types across all four de Rham spaces and the additional space `H1vec`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "spaces = {'H1': derham.V0splines, \n",
    "          'Hcurl': derham.V1splines, \n",
    "          'Hdiv': derham.V2splines, \n",
    "          'L2': derham.V3splines,\n",
    "          'H1vec': derham.Vvsplines}\n",
    "\n",
    "for name, space_info in spaces.items():\n",
    "    print(f\"\\n{name} space (spline types):\")\n",
    "    print(f\"  Spline types: {space_info.spline_types}\")\n",
    "    print(f\"  Number components: {len(space_info.quad_grid_pts)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "### 2.5 Commuting projectors\n",
    "\n",
    "The projectors that respect the de Rham sequence are accessible via the `projectors` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Commuting projectors:\")\n",
    "for space_id, proj in derham.projectors.items():\n",
    "    if space_id in derham.space_to_form.values():\n",
    "        print(f\"  P_{space_id}: {proj}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "## 3. Creating and evaluating callable spline functions\n",
    "\n",
    "Once you have a `Derham` object, you can create callable spline functions. These are wrapper objects that allow you to evaluate finite element solutions at arbitrary points in the domain."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "### 3.1 Define a test function\n",
    "\n",
    "Let's define a simple test function that we'll project into different spaces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Define a test function\n",
    "fun = lambda e1, e2, e3: np.cos(2 * np.pi * e1)\n",
    "\n",
    "# Evaluation points\n",
    "e1 = np.linspace(0, 1, 100)\n",
    "e2 = 0.5\n",
    "e3 = 0.5\n",
    "\n",
    "# Plot the function\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(e1, fun(e1, e2, e3), 'b-', linewidth=2, label='Function: cos(2π * e1)')\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('$f(e_1)$')\n",
    "plt.legend()\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "### 3.2 Project the function into the H1 space\n",
    "\n",
    "Use the `P0` projector to project the function into the discrete H1 space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project into H1 space (0-forms)\n",
    "vec = derham.P0(fun)\n",
    "print(f\"Projected vector type: {type(vec)}\")\n",
    "print(f\"Vector shape: {vec.shape}\")\n",
    "print(f\"Vector as numpy array shape: {vec[:].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "The result is a `StencilVector` (Psydac's distributed vector type). For developers, the `.toarray()` method converts it to a regular numpy array for inspection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = vec.toarray()\n",
    "print(f\"Coefficients as numpy array: shape {arr.shape}\")\n",
    "print(f\"First 10 coefficients: {arr[:10]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "### 3.3 Create a callable spline function\n",
    "\n",
    "To evaluate the finite element function at arbitrary points, create a `SplineFunction` object using `create_spline_function()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a callable spline function\n",
    "fun_h = derham.create_spline_function(\n",
    "    name=\"test_function\",\n",
    "    space_id=\"H1\",\n",
    "    coeffs=vec,\n",
    ")\n",
    "\n",
    "print(f\"Created spline function: {fun_h}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "### 3.4 Evaluate the spline function\n",
    "\n",
    "Now you can evaluate the finite element function at arbitrary points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the spline function at the same points\n",
    "fun_h_values = fun_h(e1, e2, e3, squeeze_out=True)\n",
    "\n",
    "print(\"\\nEvaluation results:\")\n",
    "print(f\"  Type: {type(fun_h_values)}\")\n",
    "print(f\"  Shape: {fun_h_values.shape}\")\n",
    "print(f\"  First 5 values: {fun_h_values[:5]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "### 3.5 Compare exact and projected functions\n",
    "\n",
    "Let's visualize the quality of the projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute exact values and projection\n",
    "exact_vals = fun(e1, e2, e3)\n",
    "proj_vals = fun_h(e1, e2, e3, squeeze_out=True)\n",
    "\n",
    "# Plot comparison\n",
    "plt.figure(figsize=(12, 5))\n",
    "\n",
    "# Left plot: comparison\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(e1, exact_vals, 'b-', linewidth=2, label='Exact')\n",
    "plt.plot(e1, proj_vals, 'r--', linewidth=2, label='H1 projection')\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('$f(e_1)$')\n",
    "plt.legend()\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Exact vs. Projected Function')\n",
    "\n",
    "# Right plot: error\n",
    "plt.subplot(1, 2, 2)\n",
    "error = np.abs(exact_vals - proj_vals)\n",
    "plt.semilogy(e1, error, 'g-', linewidth=2)\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('Absolute error')\n",
    "plt.grid(True, alpha=0.3, which='both')\n",
    "plt.title('Projection Error')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"Maximum error: {error.max():.2e}\")\n",
    "print(f\"Relative error: {error.max() / np.abs(exact_vals).max():.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "## 4. Differential operators: Grad, Curl, and Div\n",
    "\n",
    "The discrete de Rham complex contains differential operators connecting the different spaces. These are accessible as attributes of the `Derham` object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37",
   "metadata": {},
   "source": [
    "### 4.1 The Grad operator ($\\nabla$)\n",
    "\n",
    "The gradient operator maps from H1 (0-forms) to Hcurl (1-forms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access the gradient operator\n",
    "grad_op = derham.grad\n",
    "print(f\"Gradient operator: {grad_op}\")\n",
    "print(f\"Domain (input space): {grad_op.domain}\")\n",
    "print(f\"Codomain (output space): {grad_op.codomain}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39",
   "metadata": {},
   "source": [
    "Apply the gradient operator to a vector of coefficients using the `.dot()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First project a function into H1\n",
    "scalar_fun = lambda e1, e2, e3: np.cos(2 * np.pi * e1)\n",
    "scalar_vec = derham.P0(scalar_fun)\n",
    "\n",
    "# Apply gradient: H1 -> Hcurl\n",
    "grad_vec = derham.grad.dot(scalar_vec)\n",
    "print(\"\\nGradient of scalar function:\")\n",
    "print(f\"  Input type: {type(scalar_vec)}, space dimension: {scalar_vec.shape}\")\n",
    "print(f\"  Output type: {type(grad_vec)}, space dimension: {grad_vec.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41",
   "metadata": {},
   "source": [
    "Create a callable 1-form function to visualize the gradient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create callable gradient function\n",
    "grad_fun = derham.create_spline_function(\n",
    "    name=\"gradient\",\n",
    "    space_id=\"Hcurl\",\n",
    "    coeffs=grad_vec,\n",
    ")\n",
    "\n",
    "# Evaluate the gradient\n",
    "grad_vals = grad_fun(e1, e2, e3, squeeze_out=True)\n",
    "\n",
    "# The analytical gradient is -2π * sin(2π * e1)\n",
    "exact_grad = -2 * np.pi * np.sin(2 * np.pi * e1)\n",
    "\n",
    "# Plot the gradient\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(e1, exact_grad, 'b-', linewidth=2, label='Exact gradient')\n",
    "plt.plot(e1, grad_vals[0], 'r--', linewidth=2, label='Discrete gradient (e1-component)')\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('$\\\\nabla f$ value')\n",
    "plt.legend()\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Discrete Gradient')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "error_grad = np.abs(exact_grad - grad_vals[0])\n",
    "plt.semilogy(e1, error_grad, 'g-', linewidth=2)\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('Absolute error')\n",
    "plt.grid(True, alpha=0.3, which='both')\n",
    "plt.title('Gradient Error')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"Maximum gradient error: {error_grad.max():.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43",
   "metadata": {},
   "source": [
    "### 4.2 The Curl operator ($\\nabla \\times$)\n",
    "\n",
    "The curl operator maps from Hcurl (1-forms) to Hdiv (2-forms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access the curl operator\n",
    "curl_op = derham.curl\n",
    "print(f\"Curl operator: {curl_op}\")\n",
    "print(f\"Domain (input space): {curl_op.domain}\")\n",
    "print(f\"Codomain (output space): {curl_op.codomain}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45",
   "metadata": {},
   "source": [
    "To demonstrate curl, we need a vector function. Let's create one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a vector field (3 components, one for each e_i direction)\n",
    "u1 = lambda e1, e2, e3: np.sin(2 * np.pi * e1) * np.cos(2 * np.pi * e2)      # varies in e1, e2\n",
    "u2 = lambda e1, e2, e3: np.cos(2 * np.pi * e1) * np.sin(2 * np.pi * e3)      # varies in e1, e3\n",
    "u3 = lambda e1, e2, e3: np.sin(2 * np.pi * e2) * np.cos(2 * np.pi * e3)      # varies in e2, e3\n",
    "\n",
    "# Project vector field into Hcurl space\n",
    "vector_vec = derham.P1((u1, u2, u3))\n",
    "\n",
    "# Apply curl: Hcurl -> Hdiv\n",
    "curl_vec = derham.curl.dot(vector_vec)\n",
    "print(\"\\nCurl of vector field:\")\n",
    "print(f\"  Input type: {type(vector_vec)}, space dimension: {vector_vec.shape}\")\n",
    "print(f\"  Output type: {type(curl_vec)}, space dimension: {curl_vec.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47",
   "metadata": {},
   "source": [
    "### 4.3 The Div operator ($\\nabla \\cdot$)\n",
    "\n",
    "The divergence operator maps from Hdiv (2-forms) to L2 (3-forms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access the divergence operator\n",
    "div_op = derham.div\n",
    "print(f\"Divergence operator: {div_op}\")\n",
    "print(f\"Domain (input space): {div_op.domain}\")\n",
    "print(f\"Codomain (output space): {div_op.codomain}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49",
   "metadata": {},
   "source": [
    "To demonstrate divergence, we need a 2-form (which represents a vector normal to surfaces). Let's create one by projecting a vector field into Hdiv:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a vector field to project into Hdiv\n",
    "v1 = lambda e1, e2, e3: np.sin(2 * np.pi * e1) * (1 - e1)              # varies with e1\n",
    "v2 = lambda e1, e2, e3: np.cos(2 * np.pi * e2) * np.sin(2 * np.pi * e3)    # varies with e2, e3\n",
    "v3 = lambda e1, e2, e3: (1 - e1) * np.cos(2 * np.pi * e3)              # varies with e1, e3\n",
    "\n",
    "# Project into Hdiv space (2-forms)\n",
    "hdiv_vec = derham.P2((v1, v2, v3))\n",
    "\n",
    "# Apply divergence: Hdiv -> L2\n",
    "div_vec = derham.div.dot(hdiv_vec)\n",
    "print(\"\\nDivergence of 2-form:\")\n",
    "print(f\"  Input type: {type(hdiv_vec)}, space dimension: {hdiv_vec.shape}\")\n",
    "print(f\"  Output type: {type(div_vec)}, space dimension: {div_vec.shape}\")\n",
    "print(f\"  Output is a 3-form (scalar): dimension {div_vec.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51",
   "metadata": {},
   "source": [
    "## 5. Commuting projectors: P0, P1, P2, P3\n",
    "\n",
    "The geometric projectors P0, P1, P2, P3 project functions into each of the four de Rham spaces. These projectors **commute** with the differential operators, which is a fundamental property of FEEC.\n",
    "\n",
    "Let's demonstrate this with a simple 1D example:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52",
   "metadata": {},
   "source": [
    "### 5.1 Projecting a scalar function\n",
    "\n",
    "Use `P0` to project into the H1 space (0-forms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define scalar function\n",
    "scalar = lambda e1, e2, e3: 0.5 * np.sin(4 * np.pi * e1)\n",
    "\n",
    "# Project into H1\n",
    "phi_h = derham.P0(scalar)\n",
    "print(f\"P0 (H1 space):    dimension = {phi_h.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54",
   "metadata": {},
   "source": [
    "### 5.2 Projecting the gradient\n",
    "\n",
    "The analytical gradient can be projected into Hcurl using `P1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the analytical gradient\n",
    "dx = lambda e1, e2, e3: 4 * np.pi * 0.5 * np.cos(4 * np.pi * e1)  # derivative in e1 direction\n",
    "dy = lambda e1, e2, e3: np.zeros_like(e1)                                  # zero in other directions\n",
    "dz = lambda e1, e2, e3: np.zeros_like(e1)\n",
    "\n",
    "# Project gradient into Hcurl\n",
    "grad_phi_h = derham.P1((dx, dy, dz))\n",
    "print(f\"P1 (Hcurl space): dimension = {grad_phi_h.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56",
   "metadata": {},
   "source": [
    "### 5.3 Using all five projectors\n",
    "\n",
    "Here's a summary of the four projectors for the 0, 1, 2, 3-form and vector spaces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define test functions for each space\n",
    "f0 = lambda e1, e2, e3: np.sin(2 * np.pi * e1)\n",
    "f1 = lambda e1, e2, e3: np.cos(2 * np.pi * e2)\n",
    "f2 = lambda e1, e2, e3: np.cos(2 * np.pi * e3)\n",
    "\n",
    "# Apply all projectors\n",
    "v0 = derham.P0(f0)  # -> H1\n",
    "v1 = derham.P1((f0, f1, f2))  # -> Hcurl\n",
    "v2 = derham.P2((f0, f1, f2))  # -> Hdiv\n",
    "v3 = derham.P3(f0)  # -> L2\n",
    "vv = derham.Pv((f0, f1, f2))  # -> H1 vector space\n",
    "\n",
    "print(\"Projections into all de Rham spaces:\")\n",
    "print(f\"  P0 (H1):    {v0.shape}\")\n",
    "print(f\"  P1 (Hcurl): {v1.shape}\")\n",
    "print(f\"  P2 (Hdiv):  {v2.shape}\")\n",
    "print(f\"  P3 (L2):    {v3.shape}\")\n",
    "print(f\"  Pv (H1):    {vv.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58",
   "metadata": {},
   "source": [
    "## 6. The Commuting Property\n",
    "\n",
    "One of the fundamental properties of FEEC is that the geometric projectors **commute** with the differential operators. This means:\n",
    "\n",
    "$$\\mathbb P_1(\\nabla f) = \\nabla (\\mathbb P_0(f))$$\n",
    "\n",
    "for any function $f$. Let's verify this numerically:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59",
   "metadata": {},
   "source": [
    "### 6.1 Define a test function and its gradient\n",
    "\n",
    "We'll use the sinusoidal function and compute its gradient analytically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test function: f(e1) = 0.5 * sin(2π * 2 * e1)\n",
    "def fun(e1, e2, e3):\n",
    "    return 0.5 * np.sin(4 * np.pi * e1)\n",
    "\n",
    "# Analytical gradient (with zero components in other directions)\n",
    "dx = lambda e1, e2, e3: 4 * np.pi * 0.5 * np.cos(4 * np.pi * e1)\n",
    "dy = lambda e1, e2, e3: np.zeros_like(e1)\n",
    "dz = lambda e1, e2, e3: np.zeros_like(e1)\n",
    "\n",
    "print(\"Test function: f(e1) = 0.5 * sin(4π e1)\")\n",
    "print(\"Analytical gradient: ∇f = (0.5 * 4π * cos(4π e1), 0, 0)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61",
   "metadata": {},
   "source": [
    "### 6.2 Verify commutativity both ways\n",
    "\n",
    "We compute $\\nabla(P_0(f))$ and $P_1(\\nabla f)$ and check that they're equal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method 1: Project then apply gradient\n",
    "fun_h = derham.P0(fun)                        # Project function into H1\n",
    "grad_fun_h = derham.grad.dot(fun_h)           # Apply discrete gradient\n",
    "\n",
    "# Method 2: Apply gradient then project\n",
    "dfun_h = derham.P1((dx, dy, dz))                      # Project gradient into Hcurl\n",
    "\n",
    "print(\"Method 1 (project, then differentiate):\")\n",
    "print(f\"  Intermediate: P0(f) with shape {fun_h.shape}\")\n",
    "print(f\"  Result:       grad(P0(f)) with shape {grad_fun_h.shape}\")\n",
    "print()\n",
    "print(\"Method 2 (differentiate, then project):\")\n",
    "print(f\"  Result:       P1(∇f) with shape {dfun_h.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63",
   "metadata": {},
   "source": [
    "### 6.3 Compare the two methods\n",
    "\n",
    "Check that both methods give the same result by comparing their coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to arrays for comparison\n",
    "method1 = grad_fun_h.toarray()\n",
    "method2 = dfun_h.toarray()\n",
    "\n",
    "print(\"Verification of commutativity:\")\n",
    "print(f\"  Method 1 coefficients (e1 component): {method1[:10]}\")\n",
    "print(f\"  Method 2 coefficients (e1 component): {method2[:10]}\")\n",
    "print()\n",
    "\n",
    "# Check if they're numerically close\n",
    "is_commuting = np.allclose(method1, method2)\n",
    "print(f\"✓ Commuting property verified: {is_commuting}\")\n",
    "if is_commuting:\n",
    "    max_diff = np.abs(method1 - method2).max()\n",
    "    print(f\"  Maximum coefficient difference: {max_diff:.2e}\")\n",
    "else:\n",
    "    print(f\"  Maximum coefficient difference: {np.abs(method1 - method2).max():.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65",
   "metadata": {},
   "source": [
    "### 6.4 Visualize the commuting property\n",
    "\n",
    "Let's evaluate both solutions at points and plot them to see the agreement visually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create callable functions for both methods\n",
    "grad_via_method1 = derham.create_spline_function(\"grad(P0(f))\", \"Hcurl\", grad_fun_h)\n",
    "grad_via_method2 = derham.create_spline_function(\"P1(∇f)\", \"Hcurl\", dfun_h)\n",
    "\n",
    "# Evaluation points\n",
    "e1_eval = np.linspace(0, 1, 100)\n",
    "e2_eval = 0.5\n",
    "e3_eval = 0.5\n",
    "\n",
    "# Evaluate both methods\n",
    "val1 = grad_via_method1(e1_eval, e2_eval, e3_eval, squeeze_out=True)\n",
    "val2 = grad_via_method2(e1_eval, e2_eval, e3_eval, squeeze_out=True)\n",
    "dx_exact = dx(e1_eval, e2_eval, e3_eval)\n",
    "dy_exact = dy(e1_eval, e2_eval, e3_eval)\n",
    "dz_exact = dz(e1_eval, e2_eval, e3_eval)\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(12, 5))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(e1_eval, dx_exact, 'b-', linewidth=2, label='Exact gradient')\n",
    "plt.plot(e1_eval, val1[0], 'r--', linewidth=2, label='∇(P0(f))')\n",
    "plt.plot(e1_eval, val2[0], 'g:', linewidth=2, label='P1(∇f)')\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('Gradient value')\n",
    "plt.legend()\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Commuting Property: Two Computational Paths')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "diff = np.abs(val1[0] - val2[0])\n",
    "plt.semilogy(e1_eval, diff, 'k-', linewidth=2)\n",
    "plt.xlabel('$e_1$')\n",
    "plt.ylabel('Difference')\n",
    "plt.grid(True, alpha=0.3, which='both')\n",
    "plt.title('Numerical Difference Between Methods')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"Maximum difference between methods: {diff.max():.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this tutorial, you learned:\n",
    "\n",
    "1. **Derham setup** — Create a Derham object using TensorProductGrid and DerhamOptions\n",
    "2. **Attributes and SplineSpace1D** — Inspect spaces, access 1D spline information, understand the tensor product structure\n",
    "3. **Callable spline functions** — Project functions into discrete spaces and evaluate them at arbitrary points\n",
    "4. **Differential operators** — Apply grad, curl, and div operators to discrete functions\n",
    "5. **Commuting projectors** — Use P0, P1, P2, P3 to project into the four de Rham spaces\n",
    "6. **Commuting property** — Verify the fundamental FEEC property that projectors commute with derivatives\n",
    "\n",
    "These concepts form the foundation for using FEEC in Struphy for physics simulations. For boundary condition handling and more advanced topics, see the [boundary conditions tutorial](dev_tutorial_feec_bcs.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env (3.12.3)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}