{
 "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": [
    "## CeedOperator\n",
    "\n",
    "Here we show some basic examples to illustrate the `libceed.Operator` class. In libCEED, a `libceed.Operator` defines the finite/spectral element operator associated to a `libceed.QFunction` (see [the API documentation](https://libceed.org/en/latest/libCEEDapi.html#finite-element-operator-decomposition))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we create and apply a CeedOperator for the mass matrix in 1D. By applying this operator to a vector of 1's, we compute the length of this 1D domain, similar to Ex1-Volume in the [tutorial-6-shell tutorial](./tutorial-6-shell.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libceed\n",
    "import numpy as np\n",
    "\n",
    "ceed = libceed.Ceed()\n",
    "\n",
    "num_elem = 15\n",
    "p = 5\n",
    "q = 8\n",
    "num_x = num_elem + 1\n",
    "num_u = num_elem*(p-1) + 1\n",
    "\n",
    "# Vectors\n",
    "x = ceed.Vector(num_x)\n",
    "x_array = np.zeros(num_x)\n",
    "for i in range(num_x):\n",
    "  x_array[i] = i / (num_x - 1.0)\n",
    "x.set_array(x_array, cmode=libceed.USE_POINTER)\n",
    "\n",
    "q_data = ceed.Vector(num_elem*q)\n",
    "u = ceed.Vector(num_u)\n",
    "v = ceed.Vector(num_u)\n",
    "\n",
    "# Restrictions\n",
    "indices_x = np.zeros(num_x*2, dtype=\"int32\")\n",
    "for i in range(num_x):\n",
    "  indices_x[2*i+0] = i\n",
    "  indices_x[2*i+1] = i+1\n",
    "restriction_x = ceed.ElemRestriction(num_elem, 2, 1, 1, num_x, indices_x, cmode=libceed.USE_POINTER)\n",
    "\n",
    "indices_u = np.zeros(num_elem*p, dtype=\"int32\")\n",
    "for i in range(num_elem):\n",
    "  for j in range(p):\n",
    "    indices_u[p*i+j] = i*(p-1) + j\n",
    "restriction_u = ceed.ElemRestriction(num_elem, p, 1, 1, num_u, indices_u, cmode=libceed.USE_POINTER)\n",
    "strides = np.array([1, q, q], dtype=\"int32\")\n",
    "restriction_q_data = ceed.StridedElemRestriction(num_elem, q, 1, q*num_elem, strides)\n",
    "\n",
    "# Bases\n",
    "basis_x = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS)\n",
    "basis_u = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS)\n",
    "\n",
    "# QFunctions\n",
    "qf_setup = ceed.QFunctionByName(\"Mass1DBuild\")\n",
    "qf_mass = ceed.QFunctionByName(\"MassApply\")\n",
    "\n",
    "# Setup operator\n",
    "op_setup = ceed.Operator(qf_setup)\n",
    "op_setup.set_field(\"dx\", restriction_x, basis_x, libceed.VECTOR_ACTIVE)\n",
    "op_setup.set_field(\"weights\", libceed.ELEMRESTRICTION_NONE, basis_x,\n",
    "                   libceed.VECTOR_NONE)\n",
    "op_setup.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE,\n",
    "                   libceed.VECTOR_ACTIVE)\n",
    "op_setup.check()\n",
    "print('Setup operator: ', op_setup)\n",
    "\n",
    "# Mass operator\n",
    "op_mass = ceed.Operator(qf_mass)\n",
    "op_mass.set_field(\"u\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n",
    "op_mass.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE, q_data)\n",
    "op_mass.set_field(\"v\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n",
    "op_mass.check()\n",
    "print('Mass operator: ', op_mass)\n",
    "\n",
    "# Setup\n",
    "op_setup.apply(x, q_data)\n",
    "\n",
    "# Apply mass matrix\n",
    "u.set_value(1)\n",
    "op_mass.apply(u, v)\n",
    "\n",
    "# Check\n",
    "with v.array_read() as v_array:\n",
    "  print('The length of the domain is l = %4.2f'%np.sum(v_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the next example, we create and apply a CeedOperator for the Poisson operator in 1D. By applying this operator to a vector with a linear function, we compute the 'surface area' of this 1D domain, similar to Ex2-Surface in the [tutorial-6-shell tutorial](./tutorial-6-shell.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libceed\n",
    "import numpy as np\n",
    "\n",
    "ceed = libceed.Ceed()\n",
    "\n",
    "num_elem = 15\n",
    "p = 5\n",
    "q = 8\n",
    "num_x = num_elem + 1\n",
    "num_u = num_elem*(p-1) + 1\n",
    "\n",
    "# Vectors\n",
    "x = ceed.Vector(num_x)\n",
    "x_array = np.zeros(num_x)\n",
    "for i in range(num_x):\n",
    "  x_array[i] = i / (num_x - 1.0)\n",
    "x.set_array(x_array, cmode=libceed.USE_POINTER)\n",
    "\n",
    "q_data = ceed.Vector(num_elem*q)\n",
    "u = ceed.Vector(num_u)\n",
    "v = ceed.Vector(num_u)\n",
    "\n",
    "# Restrictions\n",
    "indices_x = np.zeros(num_x*2, dtype=\"int32\")\n",
    "for i in range(num_x):\n",
    "  indices_x[2*i+0] = i\n",
    "  indices_x[2*i+1] = i+1\n",
    "restriction_x = ceed.ElemRestriction(num_elem, 2, 1, 1, num_x, indices_x, cmode=libceed.USE_POINTER)\n",
    "\n",
    "indices_u = np.zeros(num_elem*p, dtype=\"int32\")\n",
    "for i in range(num_elem):\n",
    "  for j in range(p):\n",
    "    indices_u[p*i+j] = i*(p-1) + j\n",
    "restriction_u = ceed.ElemRestriction(num_elem, p, 1, 1, num_u, indices_u, cmode=libceed.USE_POINTER)\n",
    "strides = np.array([1, q, q], dtype=\"int32\")\n",
    "restriction_q_data = ceed.StridedElemRestriction(num_elem, q, 1, q*num_elem, strides)\n",
    "\n",
    "# Bases\n",
    "basis_x = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS)\n",
    "basis_u = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS)\n",
    "\n",
    "# QFunctions\n",
    "qf_setup = ceed.QFunctionByName(\"Poisson1DBuild\")\n",
    "qf_mass = ceed.QFunctionByName(\"Poisson1DApply\")\n",
    "\n",
    "# Setup operator\n",
    "op_setup = ceed.Operator(qf_setup)\n",
    "op_setup.set_field(\"dx\", restriction_x, basis_x, libceed.VECTOR_ACTIVE)\n",
    "op_setup.set_field(\"weights\", libceed.ELEMRESTRICTION_NONE, basis_x,\n",
    "                   libceed.VECTOR_NONE)\n",
    "op_setup.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE,\n",
    "                   libceed.VECTOR_ACTIVE)\n",
    "op_setup.check()\n",
    "print('Setup operator: ', op_setup)\n",
    "\n",
    "# Poisson operator\n",
    "op_poisson = ceed.Operator(qf_mass)\n",
    "op_poisson.set_field(\"du\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n",
    "op_poisson.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE, q_data)\n",
    "op_poisson.set_field(\"dv\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n",
    "op_poisson.check()\n",
    "print('Poisson operator: ', op_poisson)\n",
    "\n",
    "# Setup\n",
    "op_setup.apply(x, q_data)\n",
    "\n",
    "# Apply Poisson operator\n",
    "with u.array_write() as u_array:\n",
    "  [points, _] = ceed.lobatto_quadrature(p)\n",
    "  for elem in range(num_elem):\n",
    "      for point in range(p):\n",
    "          u_array[elem * (p - 1) + point] = (1.0 + 2.0 * elem + points[point])/(2.0 * num_elem)\n",
    "op_poisson.apply(u, v)\n",
    "\n",
    "# Check\n",
    "with v.array_read() as v_array:\n",
    "  print('The surface area of the domain is dl = %4.2f'%np.sum(abs(v_array)))"
   ]
  }
 ],
 "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
}