{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shape Function Generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get all the necessary functions for interpolation by calling the `generate_class_functions` method of any class. It returns the shape functions as symbolic matrices and also functions for fast numerical evaluation (red the corresponding docs for further information)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note** These functions are automatically generated runtime, but they can be useful if you want to build something on your own, or for educational/publication purposes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note** Generation of shape functions and their derivatives is illustrated using 8-noded hexahedron cells, but the usage applies for all cells of the library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sigmaepsilon.mesh.cells import H8\n",
    "\n",
    "shp, dshp, shpf, shpmf, dshpf = H8.Geometry.generate_class_functions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first two of the returned items are symbolic matrices for the shape functions and their derivatives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 0.125 r s t + 0.125 r s + 0.125 r t - 0.125 r + 0.125 s t - 0.125 s - 0.125 t + 0.125\\\\0.125 r s t - 0.125 r s - 0.125 r t + 0.125 r + 0.125 s t - 0.125 s - 0.125 t + 0.125\\\\- 0.125 r s t + 0.125 r s - 0.125 r t + 0.125 r - 0.125 s t + 0.125 s - 0.125 t + 0.125\\\\0.125 r s t - 0.125 r s + 0.125 r t - 0.125 r - 0.125 s t + 0.125 s - 0.125 t + 0.125\\\\0.125 r s t + 0.125 r s - 0.125 r t - 0.125 r - 0.125 s t - 0.125 s + 0.125 t + 0.125\\\\- 0.125 r s t - 0.125 r s + 0.125 r t + 0.125 r - 0.125 s t - 0.125 s + 0.125 t + 0.125\\\\0.125 r s t + 0.125 r s + 0.125 r t + 0.125 r + 0.125 s t + 0.125 s + 0.125 t + 0.125\\\\- 0.125 r s t - 0.125 r s - 0.125 r t - 0.125 r + 0.125 s t + 0.125 s + 0.125 t + 0.125\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-0.125*r*s*t + 0.125*r*s + 0.125*r*t - 0.125*r + 0.125*s*t - 0.125*s - 0.125*t + 0.125],\n",
       "[ 0.125*r*s*t - 0.125*r*s - 0.125*r*t + 0.125*r + 0.125*s*t - 0.125*s - 0.125*t + 0.125],\n",
       "[-0.125*r*s*t + 0.125*r*s - 0.125*r*t + 0.125*r - 0.125*s*t + 0.125*s - 0.125*t + 0.125],\n",
       "[ 0.125*r*s*t - 0.125*r*s + 0.125*r*t - 0.125*r - 0.125*s*t + 0.125*s - 0.125*t + 0.125],\n",
       "[ 0.125*r*s*t + 0.125*r*s - 0.125*r*t - 0.125*r - 0.125*s*t - 0.125*s + 0.125*t + 0.125],\n",
       "[-0.125*r*s*t - 0.125*r*s + 0.125*r*t + 0.125*r - 0.125*s*t - 0.125*s + 0.125*t + 0.125],\n",
       "[ 0.125*r*s*t + 0.125*r*s + 0.125*r*t + 0.125*r + 0.125*s*t + 0.125*s + 0.125*t + 0.125],\n",
       "[-0.125*r*s*t - 0.125*r*s - 0.125*r*t - 0.125*r + 0.125*s*t + 0.125*s + 0.125*t + 0.125]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 0.125 s t + 0.125 s + 0.125 t - 0.125 & - 0.125 r t + 0.125 r + 0.125 t - 0.125 & - 0.125 r s + 0.125 r + 0.125 s - 0.125\\\\0.125 s t - 0.125 s - 0.125 t + 0.125 & 0.125 r t - 0.125 r + 0.125 t - 0.125 & 0.125 r s - 0.125 r + 0.125 s - 0.125\\\\- 0.125 s t + 0.125 s - 0.125 t + 0.125 & - 0.125 r t + 0.125 r - 0.125 t + 0.125 & - 0.125 r s - 0.125 r - 0.125 s - 0.125\\\\0.125 s t - 0.125 s + 0.125 t - 0.125 & 0.125 r t - 0.125 r - 0.125 t + 0.125 & 0.125 r s + 0.125 r - 0.125 s - 0.125\\\\0.125 s t + 0.125 s - 0.125 t - 0.125 & 0.125 r t + 0.125 r - 0.125 t - 0.125 & 0.125 r s - 0.125 r - 0.125 s + 0.125\\\\- 0.125 s t - 0.125 s + 0.125 t + 0.125 & - 0.125 r t - 0.125 r - 0.125 t - 0.125 & - 0.125 r s + 0.125 r - 0.125 s + 0.125\\\\0.125 s t + 0.125 s + 0.125 t + 0.125 & 0.125 r t + 0.125 r + 0.125 t + 0.125 & 0.125 r s + 0.125 r + 0.125 s + 0.125\\\\- 0.125 s t - 0.125 s - 0.125 t - 0.125 & - 0.125 r t - 0.125 r + 0.125 t + 0.125 & - 0.125 r s - 0.125 r + 0.125 s + 0.125\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-0.125*s*t + 0.125*s + 0.125*t - 0.125, -0.125*r*t + 0.125*r + 0.125*t - 0.125, -0.125*r*s + 0.125*r + 0.125*s - 0.125],\n",
       "[ 0.125*s*t - 0.125*s - 0.125*t + 0.125,  0.125*r*t - 0.125*r + 0.125*t - 0.125,  0.125*r*s - 0.125*r + 0.125*s - 0.125],\n",
       "[-0.125*s*t + 0.125*s - 0.125*t + 0.125, -0.125*r*t + 0.125*r - 0.125*t + 0.125, -0.125*r*s - 0.125*r - 0.125*s - 0.125],\n",
       "[ 0.125*s*t - 0.125*s + 0.125*t - 0.125,  0.125*r*t - 0.125*r - 0.125*t + 0.125,  0.125*r*s + 0.125*r - 0.125*s - 0.125],\n",
       "[ 0.125*s*t + 0.125*s - 0.125*t - 0.125,  0.125*r*t + 0.125*r - 0.125*t - 0.125,  0.125*r*s - 0.125*r - 0.125*s + 0.125],\n",
       "[-0.125*s*t - 0.125*s + 0.125*t + 0.125, -0.125*r*t - 0.125*r - 0.125*t - 0.125, -0.125*r*s + 0.125*r - 0.125*s + 0.125],\n",
       "[ 0.125*s*t + 0.125*s + 0.125*t + 0.125,  0.125*r*t + 0.125*r + 0.125*t + 0.125,  0.125*r*s + 0.125*r + 0.125*s + 0.125],\n",
       "[-0.125*s*t - 0.125*s - 0.125*t - 0.125, -0.125*r*t - 0.125*r + 0.125*t + 0.125, -0.125*r*s - 0.125*r + 0.125*s + 0.125]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dshp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last three of the returned items are functions for the numerical evaluation of shape functions, the shape function matrix and the shape function derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1., -1., -1.],\n",
       "       [ 1., -1., -1.],\n",
       "       [ 1.,  1., -1.],\n",
       "       [-1.,  1., -1.],\n",
       "       [-1., -1.,  1.],\n",
       "       [ 1., -1.,  1.],\n",
       "       [ 1.,  1.,  1.],\n",
       "       [-1.,  1.,  1.]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords = H8.Geometry.master_coordinates()\n",
    "coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinates of the master element is described by a 8x3 matrix, since the cell has 8 nodes and three spatial dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 3)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate these functions using the first two coordinates of the master element to get an idea about the shapes of the resulting arrays (see the docs for the details)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 8)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shpf(coords[:2]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 1, 8)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shpmf(coords[:2]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 8, 3)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dshpf(coords[:2]).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to evaluate shape function derivatives wrt. to the local coordinate frames of cells, you need to generate a mesh first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sigmaepsilon.mesh import PolyData, PointData, CartesianFrame\n",
    "from sigmaepsilon.mesh.grid import grid\n",
    "from sigmaepsilon.mesh.cells import H8\n",
    "\n",
    "size = Lx, Ly, Lz = 800, 600, 100\n",
    "shape = nx, ny, nz = 8, 6, 2\n",
    "xbins = np.linspace(0, Lx, nx + 1)\n",
    "ybins = np.linspace(0, Ly, ny + 1)\n",
    "zbins = np.linspace(0, Lz, nz + 1)\n",
    "bins = xbins, ybins, zbins\n",
    "coords, topo = grid(bins=bins, eshape=\"H8\")\n",
    "frame = CartesianFrame(dim=3)\n",
    "\n",
    "pd = PointData(coords=coords)\n",
    "cd = H8(topo=topo, frames=frame)\n",
    "mesh = PolyData(pd, cd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observing the shape of the topology array helps to understand the shapes of the arrays of the following blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(96, 8)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topo.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the derivatives, you need the jacobian matrices of the cells evaluated at certain points. The following block evaluates the jacobian matrices of the cells at the nodes of the cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(96, 8, 3, 3)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jac = cd.jacobian_matrix()\n",
    "jac.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, you can calculate the derivatives wrt. the local frames of the cells like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(96, 2, 8, 3)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcoords = H8.Geometry.master_coordinates()\n",
    "gdshp = H8.Geometry.shape_function_derivatives(pcoords[:2], jac=jac)\n",
    "gdshp.shape"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 ('.sigeps': venv)",
   "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.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5facf25dadae24d0f6f3d9b821e9851478f51388ee31821a60476e833f1169c6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}