{
 "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": [
    "## CeedQFunction\n",
    "\n",
    "Here we show some basic examples to illustrate the `libceed.QFunction` class. In libCEED, QFunctions represent the spatial terms of the point-wise functions describing the physics at the quadrature points (see [the API documentation](https://libceed.org/en/latest/libCEEDapi.html#api-description)). As shown in the following sketch, QFunctions (such as the one depicted, which defines the Laplacian) are point-wise functions defined at quadrature points. Hence, QFunctions are independent from element shape, resolution and order.\n",
    "\n",
    "![alt text][QFunctionSchematic]\n",
    "\n",
    "[QFunctionSchematic]: ./img/QFunctionSketch.svg \"Schematic of point-wise QFunctions, defined at quadrature points, belonging to elements that can have different shape, resolution and order.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we create and view two QFunctions (for the setup and apply, respectively, of the mass operator in 1D) from the gallery of available built-in QFunctions in libCEED"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libceed\n",
    "import numpy as np\n",
    "\n",
    "ceed = libceed.Ceed()\n",
    "\n",
    "qf_setup = ceed.QFunctionByName(\"Mass1DBuild\")\n",
    "qf_mass = ceed.QFunctionByName(\"MassApply\")\n",
    "\n",
    "print(qf_setup)\n",
    "print(qf_mass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we create and evaluate a built-in identity QFunction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qf = ceed.IdentityQFunction(1, libceed.EVAL_INTERP, libceed.EVAL_INTERP)\n",
    "\n",
    "q = 8\n",
    "\n",
    "u_array = np.zeros(q, dtype=\"float64\")\n",
    "for i in range(q):\n",
    "  u_array[i] = i*i\n",
    "\n",
    "u = ceed.Vector(q)\n",
    "u.set_array(u_array, cmode=libceed.USE_POINTER)\n",
    "v = ceed.Vector(q)\n",
    "v.set_value(0)\n",
    "\n",
    "inputs = [ u ]\n",
    "outputs = [ v ]\n",
    "qf.apply(q, inputs, outputs)\n",
    "\n",
    "print('v =', v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we create and evaluate a QFunction (for the mass operator in 1D) from the gallery of available built-in QFunctions in libCEED."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qf_setup = ceed.QFunctionByName(\"Mass1DBuild\")\n",
    "qf_mass = ceed.QFunctionByName(\"MassApply\")\n",
    "\n",
    "q = 8\n",
    "\n",
    "j_array = np.zeros(q, dtype=\"float64\")\n",
    "w_array = np.zeros(q, dtype=\"float64\")\n",
    "u_array = np.zeros(q, dtype=\"float64\")\n",
    "v_true  = np.zeros(q, dtype=\"float64\")\n",
    "for i in range(q):\n",
    "  x = 2.*i/(q-1) - 1\n",
    "  j_array[i] = 1\n",
    "  w_array[i] = 1 - x*x\n",
    "  u_array[i] = 2 + 3*x + 5*x*x\n",
    "  v_true[i]  = w_array[i] * u_array[i]\n",
    "\n",
    "j = ceed.Vector(q)\n",
    "j.set_array(j_array, cmode=libceed.USE_POINTER)\n",
    "w = ceed.Vector(q)\n",
    "w.set_array(w_array, cmode=libceed.USE_POINTER)\n",
    "u = ceed.Vector(q)\n",
    "u.set_array(u_array, cmode=libceed.USE_POINTER)\n",
    "v = ceed.Vector(q)\n",
    "v.set_value(0)\n",
    "qdata = ceed.Vector(q)\n",
    "qdata.set_value(0)\n",
    "\n",
    "inputs = [ j, w ]\n",
    "outputs = [ qdata ]\n",
    "qf_setup.apply(q, inputs, outputs)\n",
    "\n",
    "inputs = [ w, u ]\n",
    "outputs = [ v ]\n",
    "qf_mass.apply(q, inputs, outputs)\n",
    "\n",
    "print('v =', v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we create and evaluate a built-in identity QFunction 3 fields per quadrature point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fields = 3\n",
    "\n",
    "qf = ceed.IdentityQFunction(fields, libceed.EVAL_INTERP, libceed.EVAL_INTERP)\n",
    "\n",
    "q = 8\n",
    "\n",
    "u_array = np.zeros(q*fields, dtype=\"float64\")\n",
    "for i in range(q*fields):\n",
    "  u_array[i] = i*i\n",
    "\n",
    "u = ceed.Vector(q*fields)\n",
    "u.set_array(u_array, cmode=libceed.USE_POINTER)\n",
    "v = ceed.Vector(q*fields)\n",
    "v.set_value(0)\n",
    "\n",
    "inputs = [ u ]\n",
    "outputs = [ v ]\n",
    "qf.apply(q, inputs, outputs)\n",
    "\n",
    "print('v =', v)"
   ]
  }
 ],
 "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
}