1{ 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "metadata": {}, 6 "source": [ 7 "# libCEED for Python examples\n", 8 "\n", 9 "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", 10 "\n", 11 "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/)." 12 ] 13 }, 14 { 15 "cell_type": "markdown", 16 "metadata": {}, 17 "source": [ 18 "## Setting up libCEED for Python\n", 19 "\n", 20 "Install libCEED for Python by running" 21 ] 22 }, 23 { 24 "cell_type": "code", 25 "execution_count": null, 26 "metadata": {}, 27 "outputs": [], 28 "source": [ 29 "! python -m pip install libceed" 30 ] 31 }, 32 { 33 "cell_type": "markdown", 34 "metadata": {}, 35 "source": [ 36 "## CeedBasis\n", 37 "\n", 38 "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))." 39 ] 40 }, 41 { 42 "cell_type": "markdown", 43 "metadata": {}, 44 "source": [ 45 "First we declare some auxiliary functions needed in the following examples" 46 ] 47 }, 48 { 49 "cell_type": "code", 50 "execution_count": null, 51 "metadata": {}, 52 "outputs": [], 53 "source": [ 54 "%matplotlib inline\n", 55 "import numpy as np\n", 56 "import matplotlib.pyplot as plt\n", 57 "plt.style.use('ggplot')\n", 58 "\n", 59 "def eval(dim, x):\n", 60 " result, center = 1, 0.1\n", 61 " for d in range(dim):\n", 62 " result *= np.tanh(x[d] - center)\n", 63 " center += 0.1\n", 64 " return result\n", 65 "\n", 66 "def feval(x_1, x_2):\n", 67 " return x_1*x_1 + x_2*x_2 + x_1*x_2 + 1\n", 68 "\n", 69 "def dfeval(x_1, x_2):\n", 70 " return 2*x_1 + x_2" 71 ] 72 }, 73 { 74 "cell_type": "markdown", 75 "metadata": {}, 76 "source": [ 77 "## $H^1$ Lagrange bases in 1D\n", 78 "\n", 79 "The Lagrange interpolation nodes are at the Gauss-Lobatto points, so interpolation to Gauss-Lobatto quadrature points is the identity." 80 ] 81 }, 82 { 83 "cell_type": "code", 84 "execution_count": null, 85 "metadata": {}, 86 "outputs": [], 87 "source": [ 88 "import libceed\n", 89 "\n", 90 "ceed = libceed.Ceed()\n", 91 "\n", 92 "b = ceed.BasisTensorH1Lagrange(\n", 93 " dim=1, # topological dimension\n", 94 " ncomp=1, # number of components\n", 95 " P=4, # number of basis functions (nodes) per dimension\n", 96 " Q=4, # number of quadrature points per dimension\n", 97 " qmode=libceed.GAUSS_LOBATTO)\n", 98 "print(b)" 99 ] 100 }, 101 { 102 "cell_type": "markdown", 103 "metadata": {}, 104 "source": [ 105 "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." 106 ] 107 }, 108 { 109 "cell_type": "code", 110 "execution_count": null, 111 "metadata": {}, 112 "outputs": [], 113 "source": [ 114 "P = b.get_num_nodes()\n", 115 "Q_viz = 50\n", 116 "basis_viz = ceed.BasisTensorH1Lagrange(1, 1, P, Q_viz, libceed.GAUSS_LOBATTO)\n", 117 "\n", 118 "# Construct P \"elements\" with one node activated\n", 119 "I = ceed.Vector(P * P)\n", 120 "with I.array_write(P, P) as x:\n", 121 " x[...] = np.eye(P)\n", 122 "\n", 123 "basis_fns = ceed.Vector(P * Q_viz)\n", 124 "basis_viz.apply(4, libceed.EVAL_INTERP, I, basis_fns)\n", 125 "\n", 126 "qpts_viz, _ = ceed.lobatto_quadrature(Q_viz)\n", 127 "with basis_fns.array_read(Q_viz, P) as B_array:\n", 128 " plt.plot(qpts_viz, B_array)\n", 129 "\n", 130 "# Mark tho Lobatto nodes\n", 131 "nodes, _ = ceed.lobatto_quadrature(P)\n", 132 "plt.plot(nodes, 0*nodes, 'ok');" 133 ] 134 }, 135 { 136 "cell_type": "markdown", 137 "metadata": {}, 138 "source": [ 139 "In contrast, the Gauss quadrature points are not collocated, and thus all basis functions are generally nonzero at every quadrature point." 140 ] 141 }, 142 { 143 "cell_type": "code", 144 "execution_count": null, 145 "metadata": {}, 146 "outputs": [], 147 "source": [ 148 "b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS)\n", 149 "print(b)\n", 150 "\n", 151 "with basis_fns.array_read(Q_viz, P) as B_array:\n", 152 " plt.plot(qpts_viz, B_array)\n", 153 "# Mark tho Gauss quadrature points\n", 154 "qpts, _ = ceed.gauss_quadrature(P)\n", 155 "plt.plot(qpts, 0*qpts, 'ok');" 156 ] 157 }, 158 { 159 "cell_type": "markdown", 160 "metadata": {}, 161 "source": [ 162 "Although the underlying functions are not an intrinsic property of a `libceed.Basis` in libCEED, the sizes are.\n", 163 "Here, we create a 3D tensor product element with more quadrature points than Lagrange interpolation nodes." 164 ] 165 }, 166 { 167 "cell_type": "code", 168 "execution_count": null, 169 "metadata": {}, 170 "outputs": [], 171 "source": [ 172 "b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO)\n", 173 "\n", 174 "p = b.get_num_nodes()\n", 175 "print('p =', p)\n", 176 "\n", 177 "q = b.get_num_quadrature_points()\n", 178 "print('q =', q)" 179 ] 180 }, 181 { 182 "cell_type": "markdown", 183 "metadata": {}, 184 "source": [ 185 "* In the following example, we demonstrate the application of an interpolatory basis in multiple dimensions" 186 ] 187 }, 188 { 189 "cell_type": "code", 190 "execution_count": null, 191 "metadata": {}, 192 "outputs": [], 193 "source": [ 194 "for dim in range(1, 4):\n", 195 " Q = 4\n", 196 " Q_dim = Q**dim\n", 197 " X_dim = 2**dim\n", 198 " x = np.empty(X_dim*dim, dtype=\"float64\")\n", 199 " u_array = np.empty(Q_dim, dtype=\"float64\")\n", 200 "\n", 201 " for d in range(dim):\n", 202 " for i in range(X_dim):\n", 203 " x[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n", 204 "\n", 205 " X = ceed.Vector(X_dim*dim)\n", 206 " X.set_array(x, cmode=libceed.USE_POINTER)\n", 207 " X_q = ceed.Vector(Q_dim*dim)\n", 208 " X_q.set_value(0)\n", 209 " U = ceed.Vector(Q_dim)\n", 210 " U.set_value(0)\n", 211 " U_q = ceed.Vector(Q_dim)\n", 212 "\n", 213 " basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)\n", 214 " basis_u_lobatto = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)\n", 215 "\n", 216 " basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n", 217 "\n", 218 " with X_q.array_read() as x_array:\n", 219 " for i in range(Q_dim):\n", 220 " x = np.empty(dim, dtype=\"float64\")\n", 221 " for d in range(dim):\n", 222 " x[d] = x_array[d*Q_dim + i]\n", 223 " u_array[i] = eval(dim, x)\n", 224 "\n", 225 " U_q.set_array(u_array, cmode=libceed.USE_POINTER)\n", 226 "\n", 227 " # This operation is the identity because the quadrature is collocated\n", 228 " basis_u_lobatto.T.apply(1, libceed.EVAL_INTERP, U_q, U)\n", 229 "\n", 230 " basis_x_gauss = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)\n", 231 " basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)\n", 232 "\n", 233 " basis_x_gauss.apply(1, libceed.EVAL_INTERP, X, X_q)\n", 234 " basis_u_gauss.apply(1, libceed.EVAL_INTERP, U, U_q)\n", 235 "\n", 236 " with X_q.array_read() as x_array, U_q.array_read() as u_array:\n", 237 " if dim == 2:\n", 238 " # Default ordering is contiguous in x direction, but\n", 239 " # pyplot expects meshgrid convention, which is transposed.\n", 240 " x, y = x_array.reshape(2, Q, Q).transpose(0, 2, 1)\n", 241 " plt.scatter(x, y, c=np.array(u_array).reshape(Q, Q))\n", 242 " plt.xlim(-1, 1)\n", 243 " plt.ylim(-1, 1)\n", 244 " plt.colorbar(label='u')" 245 ] 246 }, 247 { 248 "cell_type": "markdown", 249 "metadata": {}, 250 "source": [ 251 "* In the following example, we demonstrate the application of the gradient of the shape functions in multiple dimensions" 252 ] 253 }, 254 { 255 "cell_type": "code", 256 "execution_count": null, 257 "metadata": {}, 258 "outputs": [], 259 "source": [ 260 "for dim in range (1, 4):\n", 261 " P, Q = 8, 10\n", 262 " P_dim = P**dim\n", 263 " Q_dim = Q**dim\n", 264 " X_dim = 2**dim\n", 265 " sum_1 = sum_2 = 0\n", 266 " x_array = np.empty(X_dim*dim, dtype=\"float64\")\n", 267 " u_array = np.empty(P_dim, dtype=\"float64\")\n", 268 "\n", 269 " for d in range(dim):\n", 270 " for i in range(X_dim):\n", 271 " x_array[d*X_dim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n", 272 "\n", 273 " X = ceed.Vector(X_dim*dim)\n", 274 " X.set_array(x_array, cmode=libceed.USE_POINTER)\n", 275 " X_q = ceed.Vector(P_dim*dim)\n", 276 " X_q.set_value(0)\n", 277 " U = ceed.Vector(P_dim)\n", 278 " U_q = ceed.Vector(Q_dim*dim)\n", 279 " U_q.set_value(0)\n", 280 " Ones = ceed.Vector(Q_dim*dim)\n", 281 " Ones.set_value(1)\n", 282 " G_transpose_ones = ceed.Vector(P_dim)\n", 283 " G_transpose_ones.set_value(0)\n", 284 "\n", 285 " # Get function values at quadrature points\n", 286 " basis_x_lobatto = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)\n", 287 " basis_x_lobatto.apply(1, libceed.EVAL_INTERP, X, X_q)\n", 288 "\n", 289 " with X_q.array_read() as x_array:\n", 290 " for i in range(P_dim):\n", 291 " x = np.empty(dim, dtype=\"float64\")\n", 292 " for d in range(dim):\n", 293 " x[d] = x_array[d*P_dim + i]\n", 294 " u_array[i] = eval(dim, x)\n", 295 "\n", 296 " U.set_array(u_array, cmode=libceed.USE_POINTER)\n", 297 "\n", 298 " # Calculate G u at quadrature points, G' * 1 at dofs\n", 299 " basis_u_gauss = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)\n", 300 " basis_u_gauss.apply(1, libceed.EVAL_GRAD, U, U_q)\n", 301 " basis_u_gauss.T.apply(1, libceed.EVAL_GRAD, Ones, G_transpose_ones)\n", 302 "\n", 303 " # Check if 1' * G * u = u' * (G' * 1)\n", 304 " with G_transpose_ones.array_read() as g_array, U_q.array_read() as uq_array:\n", 305 " for i in range(P_dim):\n", 306 " sum_1 += g_array[i]*u_array[i]\n", 307 " for i in range(dim*Q_dim):\n", 308 " sum_2 += uq_array[i]\n", 309 "\n", 310 " # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n", 311 " print('1T * G * u - uT * (GT * 1) =', np.abs(sum_1 - sum_2))" 312 ] 313 } 314 ], 315 "metadata": { 316 "kernelspec": { 317 "display_name": "Python 3 (ipykernel)", 318 "language": "python", 319 "name": "python3" 320 }, 321 "language_info": { 322 "codemirror_mode": { 323 "name": "ipython", 324 "version": 3 325 }, 326 "file_extension": ".py", 327 "mimetype": "text/x-python", 328 "name": "python", 329 "nbconvert_exporter": "python", 330 "pygments_lexer": "ipython3", 331 "version": "3.13.2" 332 } 333 }, 334 "nbformat": 4, 335 "nbformat_minor": 4 336} 337