{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# libCEED for Python examples\n",
    "\n",
    "This is a tutorial to illustrate the main feautures of the Python interface for [libCEED](https://github.com/CEED/libCEED/), the low-level API library for efficient high-order discretization methods developed by the co-design [Center for Efficient Exascale Discretizations](https://ceed.exascaleproject.org/) (CEED) of the [Exascale Computing Project](https://www.exascaleproject.org/) (ECP).\n",
    "\n",
    "While libCEED's focus is on high-order finite/spectral element method implementations, the approach is mostly algebraic and thus applicable to other discretizations in factored form, as explained in the [user manual](https://libceed.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up libCEED for Python\n",
    "\n",
    "Install libCEED for Python by running"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! python -m pip install libceed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CeedBasis\n",
    "\n",
    "Here we show some basic examples to illustrate the `libceed.Basis` class. In libCEED, a `libceed.Basis` defines the finite element basis and associated quadrature rule (see [the API documentation](https://libceed.org/en/latest/libCEEDapi.html#finite-element-operator-decomposition))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we declare some auxiliary functions needed in the following examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "def eval(dim, x):\n",
    "  result, center = 1, 0.1\n",
    "  for d in range(dim):\n",
    "    result *= np.tanh(x[d] - center)\n",
    "    center += 0.1\n",
    "  return result\n",
    "\n",
    "def feval(x_1, x_2):\n",
    "  return x_1*x_1 + x_2*x_2 + x_1*x_2 + 1\n",
    "\n",
    "def dfeval(x_1, x_2):\n",
    "  return 2*x_1 + x_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $H^1$ Lagrange bases in 1D\n",
    "\n",
    "The Lagrange interpolation nodes are at the Gauss-Lobatto points, so interpolation to Gauss-Lobatto quadrature points is the identity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libceed\n",
    "\n",
    "ceed = libceed.Ceed()\n",
    "\n",
    "b = ceed.BasisTensorH1Lagrange(\n",
    "    dim=1,   # topological dimension\n",
    "    ncomp=1, # number of components\n",
    "    P=4,     # number of basis functions (nodes) per dimension\n",
    "    Q=4,     # number of quadrature points per dimension\n",
    "    qmode=libceed.GAUSS_LOBATTO)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although a `libceed.Basis` is fully discrete, we can use the Lagrange construction to extend the basis to continuous functions by applying `EVAL_INTERP` to the identity.  This is the Vandermonde matrix of the continuous basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = b.get_num_nodes()\n",
    "Q_viz = 50\n",
    "basis_viz = ceed.BasisTensorH1Lagrange(1, 1, P, Q_viz, libceed.GAUSS_LOBATTO)\n",
    "\n",
    "# Construct P \"elements\" with one node activated\n",
    "I = ceed.Vector(P * P)\n",
    "with I.array_write(P, P) as x:\n",
    "    x[...] = np.eye(P)\n",
    "\n",
    "basis_fns = ceed.Vector(P * Q_viz)\n",
    "basis_viz.apply(4, libceed.EVAL_INTERP, I, basis_fns)\n",
    "\n",
    "qpts_viz, _ = ceed.lobatto_quadrature(Q_viz)\n",
    "with basis_fns.array_read(Q_viz, P) as B_array:\n",
    "    plt.plot(qpts_viz, B_array)\n",
    "\n",
    "# Mark tho Lobatto nodes\n",
    "nodes, _ = ceed.lobatto_quadrature(P)\n",
    "plt.plot(nodes, 0*nodes, 'ok');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast, the Gauss quadrature points are not collocated, and thus all basis functions are generally nonzero at every quadrature point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS)\n",
    "print(b)\n",
    "\n",
    "with basis_fns.array_read(Q_viz, P) as B_array:\n",
    "    plt.plot(qpts_viz, B_array)\n",
    "# Mark tho Gauss quadrature points\n",
    "qpts, _ = ceed.gauss_quadrature(P)\n",
    "plt.plot(qpts, 0*qpts, 'ok');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the underlying functions are not an intrinsic property of a `libceed.Basis` in libCEED, the sizes are.\n",
    "Here, we create a 3D tensor product element with more quadrature points than Lagrange interpolation nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO)\n",
    "\n",
    "p = b.get_num_nodes()\n",
    "print('p =', p)\n",
    "\n",
    "q = b.get_num_quadrature_points()\n",
    "print('q =', q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we demonstrate the application of an interpolatory basis in multiple dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dim in range(1, 4):\n",
    "  Q = 4\n",
    "  Q_dim = Q**dim\n",
    "  X_dim = 2**dim\n",
    "  x = np.empty(X_dim*dim, dtype=\"float64\")\n",
    "  u_array = np.empty(Q_dim, dtype=\"float64\")\n",
    "\n",
    "  for d in range(dim):\n",
    "    for i in range(X_dim):\n",
    "      x[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
    "\n",
    "  X = ceed.Vector(X_dim*dim)\n",
    "  X.set_array(x, cmode=libceed.USE_POINTER)\n",
    "  X_q = ceed.Vector(Q_dim*dim)\n",
    "  X_q.set_value(0)\n",
    "  U = ceed.Vector(Q_dim)\n",
    "  U.set_value(0)\n",
    "  U_q = ceed.Vector(Q_dim)\n",
    "\n",
    "  basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)\n",
    "  basis_u_lobatto = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)\n",
    "\n",
    "  basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
    "\n",
    "  with X_q.array_read() as x_array:\n",
    "    for i in range(Q_dim):\n",
    "      x = np.empty(dim, dtype=\"float64\")\n",
    "      for d in range(dim):\n",
    "        x[d] = x_array[d*Q_dim + i]\n",
    "      u_array[i] = eval(dim, x)\n",
    "\n",
    "  U_q.set_array(u_array, cmode=libceed.USE_POINTER)\n",
    "\n",
    "  # This operation is the identity because the quadrature is collocated\n",
    "  basis_u_lobatto.T.apply(1, libceed.EVAL_INTERP, U_q, U)\n",
    "\n",
    "  basis_x_gauss = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)\n",
    "  basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)\n",
    "\n",
    "  basis_x_gauss.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
    "  basis_u_gauss.apply(1, libceed.EVAL_INTERP, U, U_q)\n",
    "\n",
    "  with X_q.array_read() as x_array, U_q.array_read() as u_array:\n",
    "    if dim == 2:\n",
    "        # Default ordering is contiguous in x direction, but\n",
    "        # pyplot expects meshgrid convention, which is transposed.\n",
    "        x, y = x_array.reshape(2, Q, Q).transpose(0, 2, 1)\n",
    "        plt.scatter(x, y, c=np.array(u_array).reshape(Q, Q))\n",
    "        plt.xlim(-1, 1)\n",
    "        plt.ylim(-1, 1)\n",
    "        plt.colorbar(label='u')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we demonstrate the application of the gradient of the shape functions in multiple dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dim in range (1, 4):\n",
    "  P, Q = 8, 10\n",
    "  P_dim = P**dim\n",
    "  Q_dim = Q**dim\n",
    "  X_dim = 2**dim\n",
    "  sum_1 = sum_2 = 0\n",
    "  x_array = np.empty(X_dim*dim, dtype=\"float64\")\n",
    "  u_array = np.empty(P_dim, dtype=\"float64\")\n",
    "\n",
    "  for d in range(dim):\n",
    "    for i in range(X_dim):\n",
    "      x_array[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
    "\n",
    "  X = ceed.Vector(X_dim*dim)\n",
    "  X.set_array(x_array, cmode=libceed.USE_POINTER)\n",
    "  X_q = ceed.Vector(P_dim*dim)\n",
    "  X_q.set_value(0)\n",
    "  U = ceed.Vector(P_dim)\n",
    "  U_q = ceed.Vector(Q_dim*dim)\n",
    "  U_q.set_value(0)\n",
    "  Ones = ceed.Vector(Q_dim*dim)\n",
    "  Ones.set_value(1)\n",
    "  G_transpose_ones = ceed.Vector(P_dim)\n",
    "  G_transpose_ones.set_value(0)\n",
    "\n",
    "  # Get function values at quadrature points\n",
    "  basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)\n",
    "  basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n",
    "\n",
    "  with X_q.array_read() as x_array:\n",
    "    for i in range(P_dim):\n",
    "      x = np.empty(dim, dtype=\"float64\")\n",
    "      for d in range(dim):\n",
    "        x[d] = x_array[d*P_dim + i]\n",
    "      u_array[i] = eval(dim, x)\n",
    "\n",
    "  U.set_array(u_array, cmode=libceed.USE_POINTER)\n",
    "\n",
    "  # Calculate G u at quadrature points, G' * 1 at dofs\n",
    "  basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)\n",
    "  basis_u_gauss.apply(1, libceed.EVAL_GRAD, U, U_q)\n",
    "  basis_u_gauss.T.apply(1, libceed.EVAL_GRAD, Ones, G_transpose_ones)\n",
    "\n",
    "  # Check if 1' * G * u = u' * (G' * 1)\n",
    "  with G_transpose_ones.array_read() as g_array, U_q.array_read() as uq_array:\n",
    "    for i in range(P_dim):\n",
    "      sum_1 += g_array[i]*u_array[i]\n",
    "    for i in range(dim*Q_dim):\n",
    "      sum_2 += uq_array[i]\n",
    "\n",
    "  # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n",
    "  print('1T * G * u - uT * (GT * 1) =', np.abs(sum_1 - sum_2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}