xref: /libCEED/examples/python/tutorial-3-basis.ipynb (revision 13964f0727a62e5421e6d3b433e838b96a9ce891)

{
 "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(x1, x2):\n",
    "  return x1*x1 + x2*x2 + x1*x2 + 1\n",
    "\n",
    "def dfeval(x1, x2):\n",
    "  return 2*x1 + x2"
   ]
  },
  {
   "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",
    "nviz = 50\n",
    "bviz = ceed.BasisTensorH1Lagrange(1, 1, P, nviz, libceed.GAUSS_LOBATTO)\n",
    "\n",
    "# Construct P \"elements\" with one node activated\n",
    "I = ceed.Vector(P * P)\n",
    "with I.array(P, P) as x:\n",
    "    x[...] = np.eye(P)\n",
    "\n",
    "Bvander = ceed.Vector(P * nviz)\n",
    "bviz.apply(4, libceed.EVAL_INTERP, I, Bvander)\n",
    "\n",
    "qviz, _weight = ceed.lobatto_quadrature(nviz)\n",
    "with Bvander.array_read(nviz, P) as B:\n",
    "    plt.plot(qviz, B)\n",
    "\n",
    "# Mark tho Lobatto nodes\n",
    "qb, _weight = ceed.lobatto_quadrature(P)\n",
    "plt.plot(qb, 0*qb, '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 Bvander.array_read(nviz, P) as B:\n",
    "    plt.plot(qviz, B)\n",
    "# Mark tho Gauss quadrature points\n",
    "qb, _weight = ceed.gauss_quadrature(P)\n",
    "plt.plot(qb, 0*qb, '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",
    "  Qdim = Q**dim\n",
    "  Xdim = 2**dim\n",
    "  x = np.empty(Xdim*dim, dtype=\"float64\")\n",
    "  uq = np.empty(Qdim, dtype=\"float64\")\n",
    "\n",
    "  for d in range(dim):\n",
    "    for i in range(Xdim):\n",
    "      x[d*Xdim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
    "\n",
    "  X = ceed.Vector(Xdim*dim)\n",
    "  X.set_array(x, cmode=libceed.USE_POINTER)\n",
    "  Xq = ceed.Vector(Qdim*dim)\n",
    "  Xq.set_value(0)\n",
    "  U = ceed.Vector(Qdim)\n",
    "  U.set_value(0)\n",
    "  Uq = ceed.Vector(Qdim)\n",
    "\n",
    "  bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)\n",
    "  bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)\n",
    "\n",
    "  bxl.apply(1, libceed.EVAL_INTERP, X, Xq)\n",
    "\n",
    "  with Xq.array_read() as xq:\n",
    "    for i in range(Qdim):\n",
    "      xx = np.empty(dim, dtype=\"float64\")\n",
    "      for d in range(dim):\n",
    "        xx[d] = xq[d*Qdim + i]\n",
    "      uq[i] = eval(dim, xx)\n",
    "\n",
    "  Uq.set_array(uq, cmode=libceed.USE_POINTER)\n",
    "\n",
    "  # This operation is the identity because the quadrature is collocated\n",
    "  bul.T.apply(1, libceed.EVAL_INTERP, Uq, U)\n",
    "\n",
    "  bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)\n",
    "  bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)\n",
    "\n",
    "  bxg.apply(1, libceed.EVAL_INTERP, X, Xq)\n",
    "  bug.apply(1, libceed.EVAL_INTERP, U, Uq)\n",
    "\n",
    "  with Xq.array_read() as xq, Uq.array_read() as u:\n",
    "    #print('xq =', xq)\n",
    "    #print('u =', u)\n",
    "    if dim == 2:\n",
    "        # Default ordering is contiguous in x direction, but\n",
    "        # pyplot expects meshgrid convention, which is transposed.\n",
    "        x, y = xq.reshape(2, Q, Q).transpose(0, 2, 1)\n",
    "        plt.scatter(x, y, c=np.array(u).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",
    "  Pdim = P**dim\n",
    "  Qdim = Q**dim\n",
    "  Xdim = 2**dim\n",
    "  sum1 = sum2 = 0\n",
    "  x = np.empty(Xdim*dim, dtype=\"float64\")\n",
    "  u = np.empty(Pdim, dtype=\"float64\")\n",
    "\n",
    "  for d in range(dim):\n",
    "    for i in range(Xdim):\n",
    "      x[d*Xdim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n",
    "\n",
    "  X = ceed.Vector(Xdim*dim)\n",
    "  X.set_array(x, cmode=libceed.USE_POINTER)\n",
    "  Xq = ceed.Vector(Pdim*dim)\n",
    "  Xq.set_value(0)\n",
    "  U = ceed.Vector(Pdim)\n",
    "  Uq = ceed.Vector(Qdim*dim)\n",
    "  Uq.set_value(0)\n",
    "  Ones = ceed.Vector(Qdim*dim)\n",
    "  Ones.set_value(1)\n",
    "  Gtposeones = ceed.Vector(Pdim)\n",
    "  Gtposeones.set_value(0)\n",
    "\n",
    "  # Get function values at quadrature points\n",
    "  bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)\n",
    "  bxl.apply(1, libceed.EVAL_INTERP, X, Xq)\n",
    "\n",
    "  with Xq.array_read() as xq:\n",
    "    for i in range(Pdim):\n",
    "      xx = np.empty(dim, dtype=\"float64\")\n",
    "      for d in range(dim):\n",
    "        xx[d] = xq[d*Pdim + i]\n",
    "      u[i] = eval(dim, xx)\n",
    "\n",
    "  U.set_array(u, cmode=libceed.USE_POINTER)\n",
    "\n",
    "  # Calculate G u at quadrature points, G' * 1 at dofs\n",
    "  bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)\n",
    "  bug.apply(1, libceed.EVAL_GRAD, U, Uq)\n",
    "  bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones)\n",
    "\n",
    "  # Check if 1' * G * u = u' * (G' * 1)\n",
    "  with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq:\n",
    "    for i in range(Pdim):\n",
    "      sum1 += gtposeones[i]*u[i]\n",
    "    for i in range(dim*Qdim):\n",
    "      sum2 += uq[i]\n",
    "\n",
    "  # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n",
    "  print('1T * G * u - uT * (GT * 1) =', np.abs(sum1 - sum2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advanced topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we demonstrate the QR factorization of a basis matrix.\n",
    "The representation is similar to LAPACK's [`dgeqrf`](https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html#ga3766ea903391b5cf9008132f7440ec7b), with elementary reflectors in the lower triangular block, scaled by `tau`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype=\"float64\")\n",
    "tau = np.empty(3, dtype=\"float64\")\n",
    "\n",
    "qr, tau = ceed.qr_factorization(qr, tau, 4, 3)\n",
    "\n",
    "print('qr =')\n",
    "print(qr.reshape(4, 3))\n",
    "\n",
    "print('tau =')\n",
    "print(tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we demonstrate the symmetric Schur decomposition of a basis matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n",
    "              0.0745459, 1., 0.16666509, -0.07448852,\n",
    "              -0.07448852, 0.16666509, 1., 0.0745459,\n",
    "              0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n",
    "\n",
    "lam = ceed.symmetric_schur_decomposition(A, 4)\n",
    "\n",
    "print(\"Q =\")\n",
    "for i in range(4):\n",
    "  for j in range(4):\n",
    "    if A[j+4*i] <= 1E-14 and A[j+4*i] >= -1E-14:\n",
    "       A[j+4*i] = 0\n",
    "    print(\"%12.8f\"%A[j+4*i])\n",
    "\n",
    "print(\"lambda =\")\n",
    "for i in range(4):\n",
    "  if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n",
    "    lam[i] = 0\n",
    "  print(\"%12.8f\"%lam[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the following example, we demonstrate the simultaneous diagonalization of a basis matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n",
    "              0.0745459, 1., 0.16666509, -0.07448852,\n",
    "              -0.07448852, 0.16666509, 1., 0.0745459,\n",
    "              0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n",
    "K = np.array([3.03344425, -3.41501767, 0.49824435, -0.11667092,\n",
    "              -3.41501767, 5.83354662, -2.9167733, 0.49824435,\n",
    "              0.49824435, -2.9167733, 5.83354662, -3.41501767,\n",
    "              -0.11667092, 0.49824435, -3.41501767, 3.03344425], dtype=\"float64\")\n",
    "\n",
    "x, lam = ceed.simultaneous_diagonalization(K, M, 4)\n",
    "\n",
    "print(\"x =\")\n",
    "for i in range(4):\n",
    "  for j in range(4):\n",
    "    if x[j+4*i] <= 1E-14 and x[j+4*i] >= -1E-14:\n",
    "      x[j+4*i] = 0\n",
    "    print(\"%12.8f\"%x[j+4*i])\n",
    "\n",
    "print(\"lambda =\")\n",
    "for i in range(4):\n",
    "  if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n",
    "    lam[i] = 0\n",
    "  print(\"%12.8f\"%lam[i])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}