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.readthedocs.io/)." 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 "## CeedElemRestriction\n", 37 "\n", 38 "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.readthedocs.io/en/latest/libCEEDapi.html#finite-element-operator-decomposition)).\n", 39 "\n", 40 "Here we illustrate the simple creation and application of a `libceed.ElemRestriction`, with user provided dof indices" 41 ] 42 }, 43 { 44 "cell_type": "code", 45 "execution_count": null, 46 "metadata": {}, 47 "outputs": [], 48 "source": [ 49 "import libceed\n", 50 "import numpy as np\n", 51 "\n", 52 "# In this 1D example, the dofs are indexed as\n", 53 "#\n", 54 "# Restriction input:\n", 55 "# x -- x -- x -- x\n", 56 "# 10 -- 11 -- 12 -- 13\n", 57 "#\n", 58 "# Restriction output:\n", 59 "# x -- x | x -- x | x -- x\n", 60 "# 10 -- 11 | 11 -- 12 | 12 -- 13\n", 61 "\n", 62 "ceed = libceed.Ceed()\n", 63 "\n", 64 "ne = 3\n", 65 "\n", 66 "x = ceed.Vector(ne+1)\n", 67 "a = np.arange(10, 10 + ne+1, dtype=\"float64\")\n", 68 "x.set_array(a, cmode=libceed.USE_POINTER)\n", 69 "\n", 70 "ind = np.zeros(2*ne, dtype=\"int32\")\n", 71 "for i in range(ne):\n", 72 " ind[2*i+0] = i\n", 73 " ind[2*i+1] = i+1\n", 74 " \n", 75 "r = ceed.ElemRestriction(ne, 2, 1, 1, ne+1, ind, cmode=libceed.USE_POINTER)\n", 76 "\n", 77 "y = ceed.Vector(2*ne)\n", 78 "y.set_value(0)\n", 79 "\n", 80 "r.apply(x, y)\n", 81 "\n", 82 "with y.array_read() as y_array:\n", 83 " print('y =', y_array)\n" 84 ] 85 }, 86 { 87 "cell_type": "markdown", 88 "metadata": {}, 89 "source": [ 90 "* In the following example, we illustrate how to extract the multiplicity of indices in an element restriction" 91 ] 92 }, 93 { 94 "cell_type": "code", 95 "execution_count": null, 96 "metadata": {}, 97 "outputs": [], 98 "source": [ 99 "# In this 1D example, there are four nodes per element\n", 100 "# \n", 101 "# x -- o -- o -- x -- o -- o -- x -- o -- o -- x\n", 102 "\n", 103 "ne = 3\n", 104 "\n", 105 "ind = np.zeros(4*ne, dtype=\"int32\")\n", 106 "\n", 107 "for i in range(ne):\n", 108 " ind[4*i+0] = i*3+0\n", 109 " ind[4*i+1] = i*3+1\n", 110 " ind[4*i+2] = i*3+2\n", 111 " ind[4*i+3] = i*3+3\n", 112 "\n", 113 "r = ceed.ElemRestriction(ne, 4, 1, 1, 3*ne+1, ind, cmode=libceed.USE_POINTER)\n", 114 "\n", 115 "mult = r.get_multiplicity()\n", 116 "\n", 117 "with mult.array_read() as m_array:\n", 118 " print('mult =', m_array)\n" 119 ] 120 }, 121 { 122 "cell_type": "markdown", 123 "metadata": {}, 124 "source": [ 125 "* 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.readthedocs.io/en/latest/libCEEDapi.html#finite-element-operator-decomposition) format, such as in Discontinuous Galerkin approximations." 126 ] 127 }, 128 { 129 "cell_type": "code", 130 "execution_count": null, 131 "metadata": {}, 132 "outputs": [], 133 "source": [ 134 "# In this 1D example, the dofs are indexed as\n", 135 "#\n", 136 "# Restriction input:\n", 137 "# x -- x -- x\n", 138 "# 10-11 -- 12-13 -- 14-15\n", 139 "#\n", 140 "# Restriction output:\n", 141 "# x -- x | x -- x | x -- x\n", 142 "# 10 -- 11 | 12 -- 13 | 14 -- 15\n", 143 "\n", 144 "ne = 3\n", 145 "\n", 146 "x = ceed.Vector(2*ne)\n", 147 "a = np.arange(10, 10 + 2*ne, dtype=\"float64\")\n", 148 "x.set_array(a, cmode=libceed.USE_POINTER)\n", 149 "\n", 150 "strides = np.array([1, 2, 2], dtype=\"int32\")\n", 151 "\n", 152 "r = ceed.StridedElemRestriction(ne, 2, 1, 2*ne, strides)\n", 153 "\n", 154 "y = ceed.Vector(2*ne)\n", 155 "y.set_value(0)\n", 156 "\n", 157 "r.apply(x, y)\n", 158 "\n", 159 "with y.array_read() as y_array:\n", 160 " print('y =', y_array)\n" 161 ] 162 }, 163 { 164 "cell_type": "markdown", 165 "metadata": {}, 166 "source": [ 167 "* In the following example, we illustrate the creation and view of a blocked strided (identity) element restriction" 168 ] 169 }, 170 { 171 "cell_type": "code", 172 "execution_count": null, 173 "metadata": {}, 174 "outputs": [], 175 "source": [ 176 "# In this 1D example, there are three elements (four nodes in total) \n", 177 "# \n", 178 "# x -- x -- x -- x\n", 179 "\n", 180 "ne = 3\n", 181 "\n", 182 "strides = np.array([1, 2, 2], dtype=\"int32\")\n", 183 "\n", 184 "r = ceed.BlockedStridedElemRestriction(ne, 2, 2, 1, ne+1, strides)\n", 185 "\n", 186 "print(r)" 187 ] 188 }, 189 { 190 "cell_type": "markdown", 191 "metadata": {}, 192 "source": [ 193 "### Advanced topics" 194 ] 195 }, 196 { 197 "cell_type": "markdown", 198 "metadata": {}, 199 "source": [ 200 "* 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)" 201 ] 202 }, 203 { 204 "cell_type": "code", 205 "execution_count": null, 206 "metadata": {}, 207 "outputs": [], 208 "source": [ 209 "# In this 1D example, the dofs are indexed as\n", 210 "# \n", 211 "# x -- x -- x -- x -- x -- x -- x -- x -- x\n", 212 "# 10 -- 11 -- 12 -- 13 -- 14 -- 15 -- 16 -- 17 -- 18\n", 213 "# \n", 214 "# We block elements into groups of 5:\n", 215 "# ________________________________________________________________________________________\n", 216 "# | block 0: | block 1: |\n", 217 "# | e0 | e1 | e2 | e3 | e4 | e0 | e1 | e2 |\n", 218 "# | | |\n", 219 "# | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x |\n", 220 "# | 10 -- 11 | 11 12 | 12 13 | 13 14 | 14 15 | 15 16 | 16 17 | 17 18 |\n", 221 "#\n", 222 "# Intermediate logical representation:\n", 223 "# ______________________________________________________________________________\n", 224 "# | block 0: | block 1: |\n", 225 "# | node0: node1: | node0: node1: |\n", 226 "# | 10-11-12-13-14 11-12-13-14-15 | 15-16-17- *- * 16-17-18- *- * |\n", 227 "# | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 |\n", 228 "#\n", 229 "# Forward restriction output:\n", 230 "# ______________________________________________________________________________\n", 231 "# | block 0: | block 1: |\n", 232 "# | node0: node1: | node0: node1: |\n", 233 "# | 10-11-12-13-14 11-12-13-14-15 | 15-16-17-17-17 16-17-18-18-18 |\n", 234 "# | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 |\n", 235 "\n", 236 "ne = 8\n", 237 "blksize = 5\n", 238 "\n", 239 "x = ceed.Vector(ne+1)\n", 240 "a = np.arange(10, 10 + ne+1, dtype=\"float64\")\n", 241 "x.set_array(a, cmode=libceed.USE_POINTER)\n", 242 "\n", 243 "ind = np.zeros(2*ne, dtype=\"int32\")\n", 244 "for i in range(ne):\n", 245 " ind[2*i+0] = i\n", 246 " ind[2*i+1] = i+1\n", 247 "\n", 248 "r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne+1, ind,\n", 249 " cmode=libceed.USE_POINTER)\n", 250 "\n", 251 "y = ceed.Vector(2*blksize*2)\n", 252 "y.set_value(0)\n", 253 "\n", 254 "r.apply(x, y)\n", 255 "\n", 256 "with y.array_read() as y_array:\n", 257 " print('y =', y_array)\n", 258 "\n", 259 "x.set_value(0)\n", 260 "r.T.apply(y, x)\n", 261 "\n", 262 "with x.array_read() as x_array:\n", 263 " print('x =', x_array)" 264 ] 265 }, 266 { 267 "cell_type": "markdown", 268 "metadata": {}, 269 "source": [ 270 "* 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)" 271 ] 272 }, 273 { 274 "cell_type": "code", 275 "execution_count": null, 276 "metadata": {}, 277 "outputs": [], 278 "source": [ 279 "# In this 1D example, the dofs are indexed as\n", 280 "# \n", 281 "# x -- x -- x -- x -- x -- x -- x -- x -- x\n", 282 "# 10 -- 11 -- 12 -- 13 -- 14 -- 15 -- 16 -- 17 -- 18\n", 283 "#\n", 284 "# We block elements into groups of 5:\n", 285 "# ________________________________________________________________________________________\n", 286 "# | block 0: | block 1: |\n", 287 "# | e0 | e1 | e2 | e3 | e4 | e0 | e1 | e2 |\n", 288 "# | | |\n", 289 "# | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x | x -- x |\n", 290 "# | 10 -- 11 | 11 12 | 12 13 | 13 14 | 14 15 | 15 16 | 16 17 | 17 18 |\n", 291 "#\n", 292 "# Intermediate logical representation (extraction of block1 only in this case):\n", 293 "# _______________________________________\n", 294 "# | block 1: |\n", 295 "# | node0: node1: |\n", 296 "# | 15-16-17- *- * 16-17-18- *- * |\n", 297 "# | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 |\n", 298 "#\n", 299 "# Forward restriction output:\n", 300 "# _______________________________________\n", 301 "# | block 1: |\n", 302 "# | node0: node1: |\n", 303 "# | 15-16-17-17-17 16-17-18-18-18 |\n", 304 "# | e0 e1 e2 e3 e4 e0 e1 e2 e3 e4 |\n", 305 "\n", 306 "ne = 8\n", 307 "blksize = 5\n", 308 "\n", 309 "x = ceed.Vector(ne+1)\n", 310 "a = np.arange(10, 10 + ne+1, dtype=\"float64\")\n", 311 "x.set_array(a, cmode=libceed.USE_POINTER)\n", 312 "\n", 313 "ind = np.zeros(2*ne, dtype=\"int32\")\n", 314 "for i in range(ne):\n", 315 " ind[2*i+0] = i\n", 316 " ind[2*i+1] = i+1\n", 317 "\n", 318 "r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne+1, ind,\n", 319 " cmode=libceed.USE_POINTER)\n", 320 "\n", 321 "y = ceed.Vector(blksize*2)\n", 322 "y.set_value(0)\n", 323 "\n", 324 "r.apply_block(1, x, y)\n", 325 "\n", 326 "with y.array_read() as y_array:\n", 327 " print('y =', y_array)\n", 328 "\n", 329 "x.set_value(0)\n", 330 "r.T.apply_block(1, y, x)\n", 331 "\n", 332 "with x.array_read() as array:\n", 333 " print('x =', x_array)" 334 ] 335 }, 336 { 337 "cell_type": "markdown", 338 "metadata": {}, 339 "source": [ 340 "Note that the nodes at the boundary between elements have multiplicty 2, while the internal nodes and the outer boundary nodes, have multiplicity 1." 341 ] 342 } 343 ], 344 "metadata": { 345 "kernelspec": { 346 "display_name": "Python 3", 347 "language": "python", 348 "name": "python3" 349 }, 350 "language_info": { 351 "codemirror_mode": { 352 "name": "ipython", 353 "version": 3 354 }, 355 "file_extension": ".py", 356 "mimetype": "text/x-python", 357 "name": "python", 358 "nbconvert_exporter": "python", 359 "pygments_lexer": "ipython3", 360 "version": "3.8.5" 361 } 362 }, 363 "nbformat": 4, 364 "nbformat_minor": 4 365} 366