{
 "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": [
    "## CeedElemRestriction\n",
    "\n",
    "Here we show some basic examples to illustrate the `libceed.ElemRestriction` class. In libCEED, a `libceed.ElemRestriction` groups the degrees of freedom (dofs) of the local vector according to the different elements they belong to (see [the API documentation](https://libceed.org/en/latest/libCEEDapi.html#finite-element-operator-decomposition)).\n",
    "\n",
    "Here we illustrate the simple creation and application of a `libceed.ElemRestriction`, with user provided dof indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libceed\n",
    "import numpy as np\n",
    "\n",
    "# In this 1D example, the dofs are indexed as\n",
    "#\n",
    "# Restriction input:\n",
    "#  x --  x --  x --  x\n",
    "# 10 -- 11 -- 12 -- 13\n",
    "#\n",
    "# Restriction output:\n",
    "#  x --  x |  x --  x | x --  x\n",
    "# 10 -- 11 | 11 -- 12 | 12 -- 13\n",
    "\n",
    "ceed = libceed.Ceed()\n",
    "\n",
    "num_elem = 3\n",
    "\n",
    "x = ceed.Vector(num_elem+1)\n",
    "a = np.arange(10, 10 + num_elem+1, dtype=\"float64\")\n",
    "x.set_array(a, cmode=libceed.USE_POINTER)\n",
    "\n",
    "indices = np.zeros(2*num_elem, dtype=\"int32\")\n",
    "for i in range(num_elem):\n",
    "  indices[2*i+0] = i\n",
    "  indices[2*i+1] = i+1\n",
    "    \n",
    "r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem+1, indices, cmode=libceed.USE_POINTER)\n",
    "\n",
    "y = ceed.Vector(2*num_elem)\n",
    "y.set_value(0)\n",
    "\n",
    "r.apply(x, y)\n",
    "\n",
    "with y.array_read() as y_array:\n",
    "  print('y =', y_array)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we illustrate how to extract the multiplicity of indices in an element restriction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this 1D example, there are four nodes per element\n",
    "# \n",
    "#  x -- o -- o -- x -- o -- o -- x -- o -- o -- x\n",
    "\n",
    "num_elem = 3\n",
    "\n",
    "indices = np.zeros(4*num_elem, dtype=\"int32\")\n",
    "\n",
    "for i in range(num_elem):\n",
    "  indices[4*i+0] = i*3+0\n",
    "  indices[4*i+1] = i*3+1\n",
    "  indices[4*i+2] = i*3+2\n",
    "  indices[4*i+3] = i*3+3\n",
    "\n",
    "r = ceed.ElemRestriction(num_elem, 4, 1, 1, 3*num_elem+1, indices, cmode=libceed.USE_POINTER)\n",
    "\n",
    "mult = r.get_multiplicity()\n",
    "\n",
    "with mult.array_read() as m_array:\n",
    "  print('mult =', m_array)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we illustrate the creation and use of a strided (identity) element restriction. Strided restrictions are typically used for data stored at quadrature points or for vectors stored in the [E-vector](https://libceed.org/en/latest/libCEEDapi.html#finite-element-operator-decomposition) format, such as in Discontinuous Galerkin approximations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this 1D example, the dofs are indexed as\n",
    "#\n",
    "# Restriction input:\n",
    "#   x   --   x   --   x\n",
    "# 10-11 -- 12-13 -- 14-15\n",
    "#\n",
    "# Restriction output:\n",
    "#  x --  x |  x --  x |  x --  x\n",
    "# 10 -- 11 | 12 -- 13 | 14 -- 15\n",
    "\n",
    "num_elem = 3\n",
    "\n",
    "x = ceed.Vector(2*num_elem)\n",
    "a = np.arange(10, 10 + 2*num_elem, dtype=\"float64\")\n",
    "x.set_array(a, cmode=libceed.USE_POINTER)\n",
    "\n",
    "strides = np.array([1, 2, 2], dtype=\"int32\")\n",
    "\n",
    "r = ceed.StridedElemRestriction(num_elem, 2, 1, 2*num_elem, strides)\n",
    "\n",
    "y = ceed.Vector(2*num_elem)\n",
    "y.set_value(0)\n",
    "\n",
    "r.apply(x, y)\n",
    "\n",
    "with y.array_read() as y_array:\n",
    "  print('y =', y_array)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we illustrate the creation and view of a blocked strided (identity) element restriction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this 1D example, there are three elements (four nodes in total) \n",
    "# \n",
    "#  x -- x -- x -- x\n",
    "\n",
    "num_elem = 3\n",
    "\n",
    "strides = np.array([1, 2, 2], dtype=\"int32\")\n",
    "\n",
    "r = ceed.BlockedStridedElemRestriction(num_elem, 2, 2, 1, 2*(num_elem+1), strides)\n",
    "\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advanced topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example (intended for backend developers), we illustrate the creation of a blocked element restriction (from an L-vector to an E-vector) and its transpose (inverse operation, from an E-vector to an L-vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this 1D example, the dofs are indexed as\n",
    "# \n",
    "#  x --  x --  x --  x --  x --  x --  x --  x --  x\n",
    "# 10 -- 11 -- 12 -- 13 -- 14 -- 15 -- 16 -- 17 -- 18\n",
    "# \n",
    "# We block elements into groups of 5:\n",
    "#  ________________________________________________________________________________________\n",
    "# |                   block 0:                            |         block 1:               |\n",
    "# |     e0    |    e1    |    e2    |    e3    |    e4    |    e0    |    e1    |    e2    |\n",
    "# |                                                       |                                |\n",
    "# |  x --  x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x  -- x |\n",
    "# | 10 -- 11  | 11   12  | 12   13  | 13   14  | 14   15  | 15   16  | 16   17  |  17   18 |\n",
    "#\n",
    "# Intermediate logical representation:\n",
    "#  ______________________________________________________________________________\n",
    "# |               block 0:               |               block 1:               |\n",
    "# |     node0:                 node1:    |     node0:                 node1:    |\n",
    "# | 10-11-12-13-14        11-12-13-14-15 | 15-16-17- *- *        16-17-18- *- * |\n",
    "# | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 |\n",
    "#\n",
    "# Forward restriction output:\n",
    "#  ______________________________________________________________________________\n",
    "# |               block 0:               |               block 1:               |\n",
    "# |     node0:                 node1:    |     node0:                 node1:    |\n",
    "# | 10-11-12-13-14        11-12-13-14-15 | 15-16-17-17-17        16-17-18-18-18 |\n",
    "# | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 |\n",
    "\n",
    "num_elem = 8\n",
    "block_size = 5\n",
    "\n",
    "x = ceed.Vector(num_elem+1)\n",
    "a = np.arange(10, 10 + num_elem+1, dtype=\"float64\")\n",
    "x.set_array(a, cmode=libceed.USE_POINTER)\n",
    "\n",
    "indices = np.zeros(2*num_elem, dtype=\"int32\")\n",
    "for i in range(num_elem):\n",
    "  indices[2*i+0] = i\n",
    "  indices[2*i+1] = i+1\n",
    "\n",
    "r = ceed.BlockedElemRestriction(num_elem, 2, block_size, 1, 1, num_elem+1, indices,\n",
    "                                cmode=libceed.USE_POINTER)\n",
    "\n",
    "y = ceed.Vector(2*block_size*2)\n",
    "y.set_value(0)\n",
    "\n",
    "r.apply(x, y)\n",
    "\n",
    "with y.array_read() as y_array:\n",
    "  print('y =', y_array)\n",
    "\n",
    "x.set_value(0)\n",
    "r.T.apply(y, x)\n",
    "\n",
    "with x.array_read() as x_array:\n",
    "  print('x =', x_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example (intended for backend developers), we illustrate the creation and application of a blocked element restriction (from an L-vector to an E-vector) and its transpose (inverse operation, from an E-vector to an L-vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this 1D example, the dofs are indexed as\n",
    "# \n",
    "#  x --  x --  x --  x --  x --  x --  x --  x --  x\n",
    "# 10 -- 11 -- 12 -- 13 -- 14 -- 15 -- 16 -- 17 -- 18\n",
    "#\n",
    "# We block elements into groups of 5:\n",
    "#  ________________________________________________________________________________________\n",
    "# |                   block 0:                            |         block 1:               |\n",
    "# |     e0    |    e1    |    e2    |    e3    |    e4    |    e0    |    e1    |    e2    |\n",
    "# |                                                       |                                |\n",
    "# |  x --  x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x -- x  |  x  -- x |\n",
    "# | 10 -- 11  | 11   12  | 12   13  | 13   14  | 14   15  | 15   16  | 16   17  |  17   18 |\n",
    "#\n",
    "# Intermediate logical representation (extraction of block1 only in this case):\n",
    "#  _______________________________________\n",
    "# |               block 1:               |\n",
    "# |     node0:                 node1:    |\n",
    "# | 15-16-17- *- *        16-17-18- *- * |\n",
    "# | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 |\n",
    "#\n",
    "# Forward restriction output:\n",
    "#  _______________________________________\n",
    "# |               block 1:               |\n",
    "# |     node0:                 node1:    |\n",
    "# | 15-16-17-17-17        16-17-18-18-18 |\n",
    "# | e0 e1 e2 e3 e4        e0 e1 e2 e3 e4 |\n",
    "\n",
    "num_elem = 8\n",
    "block_size = 5\n",
    "\n",
    "x = ceed.Vector(num_elem+1)\n",
    "a = np.arange(10, 10 + num_elem+1, dtype=\"float64\")\n",
    "x.set_array(a, cmode=libceed.USE_POINTER)\n",
    "\n",
    "indices = np.zeros(2*num_elem, dtype=\"int32\")\n",
    "for i in range(num_elem):\n",
    "  indices[2*i+0] = i\n",
    "  indices[2*i+1] = i+1\n",
    "\n",
    "r = ceed.BlockedElemRestriction(num_elem, 2, block_size, 1, 1, num_elem+1, indices,\n",
    "                                cmode=libceed.USE_POINTER)\n",
    "\n",
    "y = ceed.Vector(block_size*2)\n",
    "y.set_value(0)\n",
    "\n",
    "r.apply_block(1, x, y)\n",
    "\n",
    "with y.array_read() as y_array:\n",
    "  print('y =', y_array)\n",
    "\n",
    "x.set_value(0)\n",
    "r.T.apply_block(1, y, x)\n",
    "\n",
    "with x.array_read() as array:\n",
    "  print('x =', x_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the nodes at the boundary between elements have multiplicty 2, while the internal nodes and the outer boundary nodes, have multiplicity 1."
   ]
  }
 ],
 "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
}