diff --git a/README.md b/README.md index 67c36346..7ae2c037 100644 --- a/README.md +++ b/README.md @@ -16,16 +16,16 @@ If you want to contribute to this tutorial, please make a fork of the repository act -j test-nightly ``` -Any code added to the tutorial should work in parallel. - Alternatively, if you want to add a separate chapter, a Jupyter notebook can be added to a pull request, without integrating it into the tutorial. If so, the notebook will be reviewed and modified to be included in the tutorial. -Also ensure that both Python file and notebook files are updated by using jupytext, i.e. +Any code added to the tutorial should work in parallel. If any changes are made to `ipynb` files, please ensure that these changes are reflected in the corresponding `py` files by using [`jupytext`](https://jupytext.readthedocs.io/en/latest/faq.html#can-i-use-jupytext-with-jupyterhub-binder-nteract-colab-saturn-or-azure): ```bash python3 -m jupytext --sync */*.ipynb ``` +Any code added to the tutorial should work in parallel. + ## Dependencies It is adviced to use a pre-installed version of DOLFINx, for instance through conda or docker. Remaining dependencies can be installed with diff --git a/chapter1/complex_mode.ipynb b/chapter1/complex_mode.ipynb index bf63c78b..c8919636 100644 --- a/chapter1/complex_mode.ipynb +++ b/chapter1/complex_mode.ipynb @@ -73,7 +73,7 @@ "source": [ "However, as we would like to solve linear algebra problems of the form $Ax=b$, we need to be able to use matrices and vectors that support real and complex numbers. As [PETSc](https://petsc.org/release/) is one of the most popular interfaces to linear algebra packages, we need to be able to work with their matrix and vector structures.\n", "\n", - "Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choose `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode).\n", + "Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choosing `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode).\n", "\n", "We check that we are using the correct installation of PETSc by inspecting the scalar type.\n" ] @@ -162,7 +162,7 @@ "id": "9efe0968-bf32-4184-85f7-4e8cc3401cfb", "metadata": {}, "source": [ - "Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to in order to assemble vectors.\n" + "Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `ufl.derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to be able to use it to assemble vectors.\n" ] }, { diff --git a/chapter1/complex_mode.py b/chapter1/complex_mode.py index b357453a..0ef9d9fc 100644 --- a/chapter1/complex_mode.py +++ b/chapter1/complex_mode.py @@ -59,7 +59,7 @@ # However, as we would like to solve linear algebra problems of the form $Ax=b$, we need to be able to use matrices and vectors that support real and complex numbers. As [PETSc](https://petsc.org/release/) is one of the most popular interfaces to linear algebra packages, we need to be able to work with their matrix and vector structures. # -# Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choose `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode). +# Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choosing `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode). # # We check that we are using the correct installation of PETSc by inspecting the scalar type. # @@ -92,7 +92,7 @@ print(L) print(L2) -# Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to in order to assemble vectors. +# Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `ufl.derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to be able to use it to assemble vectors. # J = u_c**2 * ufl.dx diff --git a/chapter1/fundamentals_code.ipynb b/chapter1/fundamentals_code.ipynb index 809af6a6..4fb1f785 100644 --- a/chapter1/fundamentals_code.ipynb +++ b/chapter1/fundamentals_code.ipynb @@ -35,7 +35,7 @@ "\n", "Inserting $u_e$ in the original boundary problem, we find that \n", "\\begin{align}\n", - " f(x,y)= -6,\\qquad u_d(x,y)=u_e(x,y)=1+x^2+2y^2,\n", + " f(x,y)= -6,\\qquad u_D(x,y)=u_e(x,y)=1+x^2+2y^2,\n", "\\end{align}\n", "regardless of the shape of the domain as long as we prescribe \n", "$u_e$ on the boundary.\n", @@ -1541,7 +1541,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" }, "vscode": { "interpreter": { diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index cec1c568..d87deafb 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -43,7 +43,7 @@ # # Inserting $u_e$ in the original boundary problem, we find that # \begin{align} -# f(x,y)= -6,\qquad u_d(x,y)=u_e(x,y)=1+x^2+2y^2, +# f(x,y)= -6,\qquad u_D(x,y)=u_e(x,y)=1+x^2+2y^2, # \end{align} # regardless of the shape of the domain as long as we prescribe # $u_e$ on the boundary. diff --git a/chapter1/membrane.md b/chapter1/membrane.md index 1d9e7ab2..f89c6e93 100644 --- a/chapter1/membrane.md +++ b/chapter1/membrane.md @@ -8,7 +8,7 @@ In this section, we will turn our attentition to a physically more relevant prob We would like to compute the deflection $D(x,y)$ of a two-dimensional, circular membrane of radius $R$, subject to a load $p$ over the membrane. The appropriate PDE model is \begin{align} - -T \nabla^2D&=p \quad\text{in }\quad \Omega=\{(x,y)\vert x^2+y^2\leq R \}. + -T \nabla^2D&=p \quad\text{in }\quad \Omega=\{(x,y)\vert x^2+y^2\leq R^2 \}. \end{align} Here, $T$ is the tension in the membrane (constant), and $p$ is the external pressure load. The boundary of the membrane has no deflection. This implies that $D=0$ is the boundary condition. We model a localized load as a Gaussian function: \begin{align} @@ -36,7 +36,7 @@ With $D_e=\frac{AR^2}{8\pi\sigma T}$ and dropping the bars we obtain the scaled \begin{align} -\nabla^2 w = 4e^{-\beta^2(x^2+(y-R_0)^2)} \end{align} -to be solved over the unit disc with $w=0$ on the boundary. Now there are only two parameters which vary the dimensionless extent of the pressure, $\beta$, and the location of the pressure peak, $R_0\in[0,1]$. As $\beta\to 0$, the solution will approach the special case $w=1-x^2-y^2$. Given a computed scaed solution $w$, the physical deflection can be computed by +to be solved over the unit disc with $w=0$ on the boundary. Now there are only two parameters which vary the dimensionless extent of the pressure, $\beta$, and the location of the pressure peak, $R_0\in[0,1]$. As $\beta\to 0$, the solution will approach the special case $w=1-x^2-y^2$. Given a computed scaled solution $w$, the physical deflection can be computed by \begin{align} D=\frac{AR^2}{8\pi\sigma T}w. \end{align} diff --git a/chapter1/membrane_code.ipynb b/chapter1/membrane_code.ipynb index 1405ace8..0ff13c02 100644 --- a/chapter1/membrane_code.ipynb +++ b/chapter1/membrane_code.ipynb @@ -19,12 +19,12 @@ "\n", "## Creating the mesh\n", "\n", - "To create the computational geometry, we use the python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it." + "To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it." ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -42,7 +42,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -60,20 +60,9 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "1" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "gdim = 2\n", "gmsh.model.addPhysicalGroup(gdim, [membrane], 1)" @@ -89,23 +78,9 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Info : Meshing 1D...\n", - "Info : Meshing curve 1 (Ellipse)\n", - "Info : Done meshing 1D (Wall 0.00226294s, CPU 0.000809s)\n", - "Info : Meshing 2D...\n", - "Info : Meshing surface 1 (Plane, Frontal-Delaunay)\n", - "Info : Done meshing 2D (Wall 0.104361s, CPU 0.085868s)\n", - "Info : 1550 nodes 3099 elements\n" - ] - } - ], + "outputs": [], "source": [ "gmsh.option.setNumber(\"Mesh.CharacteristicLengthMin\", 0.05)\n", "gmsh.option.setNumber(\"Mesh.CharacteristicLengthMax\", 0.05)\n", @@ -124,7 +99,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -147,7 +122,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -166,7 +141,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -189,7 +164,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -213,7 +188,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -231,7 +206,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -255,7 +230,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -276,742 +251,11 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": { "scrolled": true }, - "outputs": [ - { - "data": { - "application/javascript": [ - "(function(root) {\n", - " function now() {\n", - " return new Date();\n", - " }\n", - "\n", - " var force = true;\n", - " var py_version = '3.2.1'.replace('rc', '-rc.').replace('.dev', '-dev.');\n", - " var is_dev = py_version.indexOf(\"+\") !== -1 || py_version.indexOf(\"-\") !== -1;\n", - " var reloading = false;\n", - " var Bokeh = root.Bokeh;\n", - " var bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n", - "\n", - " if (typeof (root._bokeh_timeout) === \"undefined\" || force) {\n", - " root._bokeh_timeout = Date.now() + 5000;\n", - " root._bokeh_failed_load = false;\n", - " }\n", - "\n", - " function run_callbacks() {\n", - " try {\n", - " root._bokeh_onload_callbacks.forEach(function(callback) {\n", - " if (callback != null)\n", - " callback();\n", - " });\n", - " } finally {\n", - " delete root._bokeh_onload_callbacks;\n", - " }\n", - " console.debug(\"Bokeh: all callbacks have finished\");\n", - " }\n", - "\n", - " function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n", - " if (css_urls == null) css_urls = [];\n", - " if (js_urls == null) js_urls = [];\n", - " if (js_modules == null) js_modules = [];\n", - " if (js_exports == null) js_exports = {};\n", - "\n", - " root._bokeh_onload_callbacks.push(callback);\n", - "\n", - " if (root._bokeh_is_loading > 0) {\n", - " console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n", - " return null;\n", - " }\n", - " if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n", - " run_callbacks();\n", - " return null;\n", - " }\n", - " if (!reloading) {\n", - " console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n", - " }\n", - "\n", - " function on_load() {\n", - " root._bokeh_is_loading--;\n", - " if (root._bokeh_is_loading === 0) {\n", - " console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n", - " run_callbacks()\n", - " }\n", - " }\n", - " window._bokeh_on_load = on_load\n", - "\n", - " function on_error() {\n", - " console.error(\"failed to load \" + url);\n", - " }\n", - "\n", - " var skip = [];\n", - " if (window.requirejs) {\n", - " window.requirejs.config({'packages': {}, 'paths': {'vtk': 'https://cdn.jsdelivr.net/npm/vtk.js@20.0.1/vtk', 'jspanel': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/jspanel', 'jspanel-modal': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/modal/jspanel.modal', 'jspanel-tooltip': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/tooltip/jspanel.tooltip', 'jspanel-hint': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/hint/jspanel.hint', 'jspanel-layout': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/layout/jspanel.layout', 'jspanel-contextmenu': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/contextmenu/jspanel.contextmenu', 'jspanel-dock': 'https://cdn.jsdelivr.net/npm/jspanel4@4.12.0/dist/extensions/dock/jspanel.dock', 'gridstack': 'https://cdn.jsdelivr.net/npm/gridstack@7.2.3/dist/gridstack-all', 'notyf': 'https://cdn.jsdelivr.net/npm/notyf@3/notyf.min'}, 'shim': {'vtk': {'exports': 'vtk'}, 'jspanel': {'exports': 'jsPanel'}, 'gridstack': {'exports': 'GridStack'}}});\n", - " require([\"vtk\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel\"], function(jsPanel) {\n", - "\twindow.jsPanel = jsPanel\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-modal\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-tooltip\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-hint\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-layout\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-contextmenu\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"jspanel-dock\"], function() {\n", - "\ton_load()\n", - " })\n", - " require([\"gridstack\"], function(GridStack) {\n", - "\twindow.GridStack = GridStack\n", - "\ton_load()\n", - " })\n", - " require([\"notyf\"], function() {\n", - "\ton_load()\n", - " })\n", - " root._bokeh_is_loading = css_urls.length + 10;\n", - " } else {\n", - " root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n", - " }\n", - "\n", - " var existing_stylesheets = []\n", - " var links = document.getElementsByTagName('link')\n", - " for (var i = 0; i < links.length; i++) {\n", - " var link = links[i]\n", - " if (link.href != null) {\n", - "\texisting_stylesheets.push(link.href)\n", - " }\n", - " }\n", - " for (var i = 0; i < css_urls.length; i++) {\n", - " var url = css_urls[i];\n", - " if (existing_stylesheets.indexOf(url) !== -1) {\n", - "\ton_load()\n", - "\tcontinue;\n", - " }\n", - " const element = document.createElement(\"link\");\n", - " element.onload = on_load;\n", - " element.onerror = on_error;\n", - " element.rel = \"stylesheet\";\n", - " element.type = \"text/css\";\n", - " element.href = url;\n", - " console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n", - " document.body.appendChild(element);\n", - " } if (((window['vtk'] !== undefined) && (!(window['vtk'] instanceof HTMLElement))) || window.requirejs) {\n", - " var urls = ['https://cdn.holoviz.org/panel/1.2.1/dist/bundled/abstractvtkplot/vtk.js@20.0.1/vtk.js'];\n", - " for (var i = 0; i < urls.length; i++) {\n", - " skip.push(urls[i])\n", - " }\n", - " } if (((window['jsPanel'] !== undefined) && (!(window['jsPanel'] instanceof HTMLElement))) || window.requirejs) {\n", - " var urls = ['https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/jspanel.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/modal/jspanel.modal.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/tooltip/jspanel.tooltip.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/hint/jspanel.hint.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/layout/jspanel.layout.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/contextmenu/jspanel.contextmenu.js', 'https://cdn.holoviz.org/panel/1.2.1/dist/bundled/floatpanel/jspanel4@4.12.0/dist/extensions/dock/jspanel.dock.js'];\n", - " for (var i = 0; i < urls.length; i++) {\n", - " skip.push(urls[i])\n", - " }\n", - " } if (((window['GridStack'] !== undefined) && (!(window['GridStack'] instanceof HTMLElement))) || window.requirejs) {\n", - " var urls = ['https://cdn.holoviz.org/panel/1.2.1/dist/bundled/gridstack/gridstack@7.2.3/dist/gridstack-all.js'];\n", - " for (var i = 0; i < urls.length; i++) {\n", - " skip.push(urls[i])\n", - " }\n", - " } if (((window['Notyf'] !== undefined) && (!(window['Notyf'] instanceof HTMLElement))) || window.requirejs) {\n", - " var urls = ['https://cdn.holoviz.org/panel/1.2.1/dist/bundled/notificationarea/notyf@3/notyf.min.js'];\n", - " for (var i = 0; i < urls.length; i++) {\n", - " skip.push(urls[i])\n", - " }\n", - " } var existing_scripts = []\n", - " var scripts = document.getElementsByTagName('script')\n", - " for (var i = 0; i < scripts.length; i++) {\n", - " var script = scripts[i]\n", - " if (script.src != null) {\n", - "\texisting_scripts.push(script.src)\n", - " }\n", - " }\n", - " for (var i = 0; i < js_urls.length; i++) {\n", - " var url = js_urls[i];\n", - " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", - "\tif (!window.requirejs) {\n", - "\t on_load();\n", - "\t}\n", - "\tcontinue;\n", - " }\n", - " var element = document.createElement('script');\n", - " element.onload = on_load;\n", - " element.onerror = on_error;\n", - " element.async = false;\n", - " element.src = url;\n", - " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", - " document.head.appendChild(element);\n", - " }\n", - " for (var i = 0; i < js_modules.length; i++) {\n", - " var url = js_modules[i];\n", - " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", - "\tif (!window.requirejs) {\n", - "\t on_load();\n", - "\t}\n", - "\tcontinue;\n", - " }\n", - " var element = document.createElement('script');\n", - " element.onload = on_load;\n", - " element.onerror = on_error;\n", - " element.async = false;\n", - " element.src = url;\n", - " element.type = \"module\";\n", - " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", - " document.head.appendChild(element);\n", - " }\n", - " for (const name in js_exports) {\n", - " var url = js_exports[name];\n", - " if (skip.indexOf(url) >= 0 || root[name] != null) {\n", - "\tif (!window.requirejs) {\n", - "\t on_load();\n", - "\t}\n", - "\tcontinue;\n", - " }\n", - " var element = document.createElement('script');\n", - " element.onerror = on_error;\n", - " element.async = false;\n", - " element.type = \"module\";\n", - " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", - " element.textContent = `\n", - " import ${name} from \"${url}\"\n", - " window.${name} = ${name}\n", - " window._bokeh_on_load()\n", - " `\n", - " document.head.appendChild(element);\n", - " }\n", - " if (!js_urls.length && !js_modules.length) {\n", - " on_load()\n", - " }\n", - " };\n", - "\n", - " function inject_raw_css(css) {\n", - " const element = document.createElement(\"style\");\n", - " element.appendChild(document.createTextNode(css));\n", - " document.body.appendChild(element);\n", - " }\n", - "\n", - " var js_urls = [\"https://cdn.holoviz.org/panel/1.2.1/dist/bundled/abstractvtkplot/vtk.js@20.0.1/vtk.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-3.2.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.2.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.2.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.2.1.min.js\", \"https://cdn.holoviz.org/panel/1.2.1/dist/panel.min.js\"];\n", - " var js_modules = [];\n", - " var js_exports = {};\n", - " var css_urls = [];\n", - " var inline_js = [ function(Bokeh) {\n", - " Bokeh.set_log_level(\"info\");\n", - " },\n", - "function(Bokeh) {} // ensure no trailing comma for IE\n", - " ];\n", - "\n", - " function run_inline_js() {\n", - " if ((root.Bokeh !== undefined) || (force === true)) {\n", - " for (var i = 0; i < inline_js.length; i++) {\n", - " inline_js[i].call(root, root.Bokeh);\n", - " }\n", - " // Cache old bokeh versions\n", - " if (Bokeh != undefined && !reloading) {\n", - "\tvar NewBokeh = root.Bokeh;\n", - "\tif (Bokeh.versions === undefined) {\n", - "\t Bokeh.versions = new Map();\n", - "\t}\n", - "\tif (NewBokeh.version !== Bokeh.version) {\n", - "\t Bokeh.versions.set(NewBokeh.version, NewBokeh)\n", - "\t}\n", - "\troot.Bokeh = Bokeh;\n", - " }} else if (Date.now() < root._bokeh_timeout) {\n", - " setTimeout(run_inline_js, 100);\n", - " } else if (!root._bokeh_failed_load) {\n", - " console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n", - " root._bokeh_failed_load = true;\n", - " }\n", - " root._bokeh_is_initializing = false\n", - " }\n", - "\n", - " function load_or_wait() {\n", - " // Implement a backoff loop that tries to ensure we do not load multiple\n", - " // versions of Bokeh and its dependencies at the same time.\n", - " // In recent versions we use the root._bokeh_is_initializing flag\n", - " // to determine whether there is an ongoing attempt to initialize\n", - " // bokeh, however for backward compatibility we also try to ensure\n", - " // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n", - " // before older versions are fully initialized.\n", - " if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n", - " root._bokeh_is_initializing = false;\n", - " root._bokeh_onload_callbacks = undefined;\n", - " console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n", - " load_or_wait();\n", - " } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n", - " setTimeout(load_or_wait, 100);\n", - " } else {\n", - " Bokeh = root.Bokeh;\n", - " bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n", - " root._bokeh_is_initializing = true\n", - " root._bokeh_onload_callbacks = []\n", - " if (!reloading && (!bokeh_loaded || is_dev)) {\n", - "\troot.Bokeh = undefined;\n", - " }\n", - " load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n", - "\tconsole.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n", - "\trun_inline_js();\n", - " });\n", - " }\n", - " }\n", - " // Give older versions of the autoload script a head-start to ensure\n", - " // they initialize before we start loading newer version.\n", - " setTimeout(load_or_wait, 100)\n", - "}(window));" - ], - "application/vnd.holoviews_load.v0+json": "" - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/javascript": [ - "\n", - "if ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n", - " window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n", - "}\n", - "\n", - "\n", - " function JupyterCommManager() {\n", - " }\n", - "\n", - " JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n", - " if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", - " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", - " comm_manager.register_target(comm_id, function(comm) {\n", - " comm.on_msg(msg_handler);\n", - " });\n", - " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", - " window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n", - " comm.onMsg = msg_handler;\n", - " });\n", - " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", - " google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n", - " var messages = comm.messages[Symbol.asyncIterator]();\n", - " function processIteratorResult(result) {\n", - " var message = result.value;\n", - " console.log(message)\n", - " var content = {data: message.data, comm_id};\n", - " var buffers = []\n", - " for (var buffer of message.buffers || []) {\n", - " buffers.push(new DataView(buffer))\n", - " }\n", - " var metadata = message.metadata || {};\n", - " var msg = {content, buffers, metadata}\n", - " msg_handler(msg);\n", - " return messages.next().then(processIteratorResult);\n", - " }\n", - " return messages.next().then(processIteratorResult);\n", - " })\n", - " }\n", - " }\n", - "\n", - " JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n", - " if (comm_id in window.PyViz.comms) {\n", - " return window.PyViz.comms[comm_id];\n", - " } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", - " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", - " var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n", - " if (msg_handler) {\n", - " comm.on_msg(msg_handler);\n", - " }\n", - " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", - " var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n", - " comm.open();\n", - " if (msg_handler) {\n", - " comm.onMsg = msg_handler;\n", - " }\n", - " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", - " var comm_promise = google.colab.kernel.comms.open(comm_id)\n", - " comm_promise.then((comm) => {\n", - " window.PyViz.comms[comm_id] = comm;\n", - " if (msg_handler) {\n", - " var messages = comm.messages[Symbol.asyncIterator]();\n", - " function processIteratorResult(result) {\n", - " var message = result.value;\n", - " var content = {data: message.data};\n", - " var metadata = message.metadata || {comm_id};\n", - " var msg = {content, metadata}\n", - " msg_handler(msg);\n", - " return messages.next().then(processIteratorResult);\n", - " }\n", - " return messages.next().then(processIteratorResult);\n", - " }\n", - " }) \n", - " var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n", - " return comm_promise.then((comm) => {\n", - " comm.send(data, metadata, buffers, disposeOnDone);\n", - " });\n", - " };\n", - " var comm = {\n", - " send: sendClosure\n", - " };\n", - " }\n", - " window.PyViz.comms[comm_id] = comm;\n", - " return comm;\n", - " }\n", - " window.PyViz.comm_manager = new JupyterCommManager();\n", - " \n", - "\n", - "\n", - "var JS_MIME_TYPE = 'application/javascript';\n", - "var HTML_MIME_TYPE = 'text/html';\n", - "var EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\n", - "var CLASS_NAME = 'output';\n", - "\n", - "/**\n", - " * Render data to the DOM node\n", - " */\n", - "function render(props, node) {\n", - " var div = document.createElement(\"div\");\n", - " var script = document.createElement(\"script\");\n", - " node.appendChild(div);\n", - " node.appendChild(script);\n", - "}\n", - "\n", - "/**\n", - " * Handle when a new output is added\n", - " */\n", - "function handle_add_output(event, handle) {\n", - " var output_area = handle.output_area;\n", - " var output = handle.output;\n", - " if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n", - " return\n", - " }\n", - " var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n", - " var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n", - " if (id !== undefined) {\n", - " var nchildren = toinsert.length;\n", - " var html_node = toinsert[nchildren-1].children[0];\n", - " html_node.innerHTML = output.data[HTML_MIME_TYPE];\n", - " var scripts = [];\n", - " var nodelist = html_node.querySelectorAll(\"script\");\n", - " for (var i in nodelist) {\n", - " if (nodelist.hasOwnProperty(i)) {\n", - " scripts.push(nodelist[i])\n", - " }\n", - " }\n", - "\n", - " scripts.forEach( function (oldScript) {\n", - " var newScript = document.createElement(\"script\");\n", - " var attrs = [];\n", - " var nodemap = oldScript.attributes;\n", - " for (var j in nodemap) {\n", - " if (nodemap.hasOwnProperty(j)) {\n", - " attrs.push(nodemap[j])\n", - " }\n", - " }\n", - " attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n", - " newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n", - " oldScript.parentNode.replaceChild(newScript, oldScript);\n", - " });\n", - " if (JS_MIME_TYPE in output.data) {\n", - " toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n", - " }\n", - " output_area._hv_plot_id = id;\n", - " if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n", - " window.PyViz.plot_index[id] = Bokeh.index[id];\n", - " } else {\n", - " window.PyViz.plot_index[id] = null;\n", - " }\n", - " } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n", - " var bk_div = document.createElement(\"div\");\n", - " bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n", - " var script_attrs = bk_div.children[0].attributes;\n", - " for (var i = 0; i < script_attrs.length; i++) {\n", - " toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n", - " }\n", - " // store reference to server id on output_area\n", - " output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n", - " }\n", - "}\n", - "\n", - "/**\n", - " * Handle when an output is cleared or removed\n", - " */\n", - "function handle_clear_output(event, handle) {\n", - " var id = handle.cell.output_area._hv_plot_id;\n", - " var server_id = handle.cell.output_area._bokeh_server_id;\n", - " if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n", - " var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n", - " if (server_id !== null) {\n", - " comm.send({event_type: 'server_delete', 'id': server_id});\n", - " return;\n", - " } else if (comm !== null) {\n", - " comm.send({event_type: 'delete', 'id': id});\n", - " }\n", - " delete PyViz.plot_index[id];\n", - " if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n", - " var doc = window.Bokeh.index[id].model.document\n", - " doc.clear();\n", - " const i = window.Bokeh.documents.indexOf(doc);\n", - " if (i > -1) {\n", - " window.Bokeh.documents.splice(i, 1);\n", - " }\n", - " }\n", - "}\n", - "\n", - "/**\n", - " * Handle kernel restart event\n", - " */\n", - "function handle_kernel_cleanup(event, handle) {\n", - " delete PyViz.comms[\"hv-extension-comm\"];\n", - " window.PyViz.plot_index = {}\n", - "}\n", - "\n", - "/**\n", - " * Handle update_display_data messages\n", - " */\n", - "function handle_update_output(event, handle) {\n", - " handle_clear_output(event, {cell: {output_area: handle.output_area}})\n", - " handle_add_output(event, handle)\n", - "}\n", - "\n", - "function register_renderer(events, OutputArea) {\n", - " function append_mime(data, metadata, element) {\n", - " // create a DOM node to render to\n", - " var toinsert = this.create_output_subarea(\n", - " metadata,\n", - " CLASS_NAME,\n", - " EXEC_MIME_TYPE\n", - " );\n", - " this.keyboard_manager.register_events(toinsert);\n", - " // Render to node\n", - " var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n", - " render(props, toinsert[0]);\n", - " element.append(toinsert);\n", - " return toinsert\n", - " }\n", - "\n", - " events.on('output_added.OutputArea', handle_add_output);\n", - " events.on('output_updated.OutputArea', handle_update_output);\n", - " events.on('clear_output.CodeCell', handle_clear_output);\n", - " events.on('delete.Cell', handle_clear_output);\n", - " events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n", - "\n", - " OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n", - " safe: true,\n", - " index: 0\n", - " });\n", - "}\n", - "\n", - "if (window.Jupyter !== undefined) {\n", - " try {\n", - " var events = require('base/js/events');\n", - " var OutputArea = require('notebook/js/outputarea').OutputArea;\n", - " if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n", - " register_renderer(events, OutputArea);\n", - " }\n", - " } catch(err) {\n", - " }\n", - "}\n" - ], - "application/vnd.holoviews_load.v0+json": "" - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/vnd.holoviews_exec.v0+json": "", - "text/html": [ - "
\n", - "
\n", - "
\n", - "" - ] - }, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "0b886c79-d820-4291-af03-076387355839" - } - }, - "output_type": "display_data" - }, - { - "data": {}, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/vnd.holoviews_exec.v0+json": "", - "text/html": [ - "
\n", - "
\n", - "
\n", - "" - ], - "text/plain": [ - "VTKRenderWindowSynchronized(vtkXOpenGLRenderWindow, color_mappers=[LinearColorMapper(id='00a...], sizing_mode='stretch_width')" - ] - }, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "91004001-5bc1-41b1-aaac-4c2d797bc630" - } - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "from dolfinx.plot import vtk_mesh\n", "import pyvista\n", @@ -1043,97 +287,9 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": {}, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "application/vnd.holoviews_exec.v0+json": "", - "text/html": [ - "
\n", - "
\n", - "
\n", - "" - ], - "text/plain": [ - "VTKRenderWindowSynchronized(vtkXOpenGLRenderWindow, color_mappers=[LinearColorMapper(id='c7a...], sizing_mode='stretch_width')" - ] - }, - "metadata": { - "application/vnd.holoviews_exec.v0+json": { - "id": "c3e924b8-6488-4327-a1f1-ab80c8b955a9" - } - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "load_plotter = pyvista.Plotter()\n", "p_grid = pyvista.UnstructuredGrid(*vtk_mesh(Q))\n", @@ -1160,7 +316,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1184,7 +340,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1198,8 +354,8 @@ "metadata": {}, "source": [ "Now we can compute which cells the bounding box tree collides with using `dolfinx.geometry.compute_collisions_points`. This function returns a list of cells whose bounding box collide for each input point. As different points might have different number of cells, the data is stored in `dolfinx.cpp.graph.AdjacencyList_int32`, where one can access the cells for the `i`th point by calling `links(i)`.\n", - " However, as the bounding box of a cell spans more of $\\mathbb{R}^n$ than the actual cell, we check that the actual cell collides with cell \n", - " using `dolfinx.geometry.select_colliding_cells`, who measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries).\n", + "However, as the bounding box of a cell spans more of $\\mathbb{R}^n$ than the actual cell, we check that the actual cell collides with cell\n", + "using `dolfinx.geometry.select_colliding_cells`, which measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries).\n", "This function also returns an adjacency-list, as the point might align with a facet, edge or vertex that is shared between multiple cells in the mesh.\n", "\n", "Finally, we would like the code below to run in parallel, when the mesh is distributed over multiple processors. In that case, it is not guaranteed that every point in `points` is on each processor. Therefore we create a subset `points_on_proc` only containing the points found on the current processor." @@ -1207,7 +363,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1228,12 +384,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now got a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points." + "We now have a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points.\n" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1252,20 +408,9 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB40UlEQVR4nO3dd3hT5RfA8W+6KXSwKVAqCrJH2WWLbERwoIyfZSsIihZBkCFDLCiy9yyCiBNcrIoCAkU2sgSKbGjZndCmzf39ERoaupI2q5fzeZ48JDd3nNPbNIf3vvd9NYqiKAghhBBC2ImTvQMQQgghxJNNihEhhBBC2JUUI0IIIYSwKylGhBBCCGFXUowIIYQQwq6kGBFCCCGEXUkxIoQQQgi7crF3AKbQ6XRcu3YNLy8vNBqNvcMRQgghhAkURSEuLo7SpUvj5JR1+0e+KEauXbuGv7+/vcMQQgghRC5cvnyZsmXLZvl+vihGvLy8AH0y3t7eFtuvVqtl69attG3bFldXV4vt11GoPT9Qf45qzw/Un6Pa8wP15yj55V5sbCz+/v6G7/Gs5ItiJO3SjLe3t8WLEU9PT7y9vVX7C6bm/ED9Oao9P1B/jmrPD9Sfo+SXdzl1sZAOrEIIIYSwKylGhBBCCGFXUowIIYQQwq7yRZ8RIYQQxlJTU9FqtfYOA9D3OXBxceHBgwekpqbaOxyLk/yy5urqirOzc55jkGJECCHyEUVRiIqK4t69e/YOxUBRFEqVKsXly5dVORaU5Jc9X19fSpUqlaefjRQjQgiRj6QVIiVKlMDT09Mhvhx1Oh3x8fEUKlQo24Gt8ivJL3OKopCYmMiNGzcA8PPzy3UMUowIIUQ+kZqaaihEihYtau9wDHQ6HcnJyXh4eKj2y1ryy1yBAgUAuHHjBiVKlMj1JRv1/VSFEEKl0vqIeHp62jkSIR5J+33MSx8mKUaEECKfcYRLM0KkscTvoxQjQgghhLCrPBUjU6dORaPR8N5772W73nfffUflypXx8PCgRo0abNy4MS+HFUIIIYSK5LoY2b9/P4sXL6ZmzZrZrrdnzx569OhB//79OXz4MF27dqVr164cP348t4cWQgghhIrkqhiJj4+nV69eLF26lMKFC2e77uzZs2nfvj0jRoygSpUqTJ48mTp16jBv3rxcBSyEEGoVGQlhYRATY+9IhLCtXN3aO2TIEDp16kTr1q355JNPsl03IiKCkJAQo2Xt2rVjw4YNWW6TlJREUlKS4XVsbCyg76lryREH0/blKKMYWpra8wP156j2/ED9OZqaX0oKtGrlwuXLGrZv17F0acaRMLVaLYqioNPp0Ol0Vok3NxRFMfxri7gURWHQoEH88MMP3L17l4MHDxISEkKtWrWYOXOmVY4H8Nxzz1G7dm2rHONxt2/fplq1auzdu5ennnrKqsfK6vz16NGD+vXrZ/gOf5xOp0NRFLRabYZbe039XJtdjKxbt45Dhw6xf/9+k9aPioqiZMmSRstKlixJVFRUltuEhoYyceLEDMu3bt1qlVvawsPDLb5PR6L2/ED9Oao9P1B/jjnlFxnpy+XLLQBYtcqJLl1+4fEhH1xcXChVqhTx8fEkJydbK9Rci4uLy/K9t99+m6+//hrQ51G4cGGqVavGK6+8Qs+ePc0a3yI8PJxVq1bxyy+/8NRTT1G0aFFSUlJITk42/Oc1r1544QVq1KhBaGioYVlYWBguLi4WO0Z2JkyYQIcOHShSpEiujjd16lSmTZtmtKxixYrs27fP8Hrp0qXMnTuXGzduUL16daZNm0bdunUN7w8bNoxOnTrRrVs3fHx8sjxWcnIy9+/fZ+fOnaSkpBi9l5iYaFK8ZhUjly9fZtiwYYSHh+Ph4WHOpmYZPXq0USUWGxuLv78/bdu2xdvb22LH0Wq1hIeH06ZNG1xdXS22X0eh9vxA/TmqPT9Qf46m5peUBFOnKty6pb9Nsly5jjzeJe/BgwdcvnyZQoUKWfVvsLkURSEuLg4vL68sb/N0dXWlXbt2rFixgtTUVKKjo9myZQujR4/mt99+46effsLFxbSvpKioKPz8/GjTpo1hmYuLC25ubhb7jki/v7T8ypUrZ5PbqhMTE1mzZg2bNm0yymf37t00aNAgw+/RyZMnKVq0qNF//N3d3alWrRpbt241yiltf9988w1jx45lwYIFNGjQgC+++IJXX32VU6dOUaJECQAaNWpEhQoV+OWXX3j77bezjPfBgwcUKFCA5s2bZ/i9NLWQMqsYOXjwIDdu3KBOnTqGZampqezcuZN58+aRlJSUoYmmVKlSREdHGy2Ljo6mVKlSWR7H3d0dd3f3DMtdXV2t8sfKWvt1FGrPD9Sfo9rzA/XnmFN+rq4wZgy8/77+9Z49rqT7Tyqg/3ur0WhwcnIytCTUq1cv25ZmaylVqhQHDhwAMDTtp8WWGY1Gg4eHB6VLlwbA39+fevXqERQUxPPPP8+XX37JgAED0Ol0TJs2jSVLlhAVFcWzzz7LuHHjePXVVwHo06cPq1atAsDZ2ZmAgAAuXLhgdPyc9pEW8/Tp01myZAmXL1+mZMmSvPXWW4wZM4Y+ffqwY8cOduzYwZw5cwA4evQo7777LrVr12bWrFmAvkvBiBEjWLduHbGxsdSrV4+ZM2dSv359AFq2bEnNmjXx8PBg2bJluLm5MWjQICZMmJDtz3bz5s24u7vTuHFjo3jfeecdKlasyLp16wzftadPn6Z169aEhIQwcuRIo5+3i4uL4ef9uFmzZjFw4ED69++PTqdjxowZhIeHExYWxqhRowzrde7cmW+++YahQ4dmGa+TkxMajSbT33FTP9NmFSPPP/88x44dM1rWt29fKleuzIcffpjpMLBBQUFs27bN6Pbf8PBwgoKCzDm0EEKoXsuWj57v2AHvvJPzNlFRUVy9etVqMVlbq1atqFWrFj/++CMDBgwgNDSUNWvWsGjRIipWrMjOnTv53//+R/HixWnRogWzZ8/mmWeeYcmSJezfvz/T752c9gH6FvilS5cyc+ZMmjZtyvXr1/n3338B/Y0XZ86coXr16kyaNAmdTpfpf5BHjhzJDz/8wKpVqwgICOCzzz6jXbt2REZGUqRIEQBWrVpFSEgIf//9NxEREfTp04cmTZoYteo87q+//jK6XAL6L/yNGzfSvHlzgoODWb16NefPn6dVq1Z07drVqBBJc/bsWUqXLo2HhwdBQUGEhoZSrlw5kpOTOXjwIKNHjzba//PPP09ERITRPho0aMCUKVNISkrK9GdgKWYVI15eXlSvXt1oWcGCBSlatKhheXBwMGXKlDFcZxs2bBgtWrTgiy++oFOnTqxbt44DBw6wZMkSC6UghBDqUKMG+PrCvXv6YkRRIKerAtm1MluTJY9buXJl/vnnH5KSkvj000/5/fffDf9hffrpp9m1axeLFy+mRYsW+Pj44OXlhbOzc6YxmLKPuLg4Zs+ezbx58+jduzcAzzzzDE2bNgXAx8cHNzc3PD09KVWqFDqdLsPlhoSEBBYuXEhYWBgdOnQA9H0wwsPDWb58OSNGjACgZs2afPzxx4C+z8a8efPYtm1btsXIxYsXM23RKF26NH/88QfNmjWjZ8+eRERE0Lp1axYuXJhh3YYNGxIWFkalSpW4fv06EydOpFmzZhw/fpy4uDhSU1Mz7c95+vTpDMdMTk4mKiqKgICALGPOK4tPlHfp0iWjZrrGjRuzdu1axo4dy0cffUTFihXZsGFDhqJGCCGeVNu3w9Gj0LQpBAXBpk1w6xacPAnVqmW/bdqlkvxMURQ0Gg2RkZEkJiZm+KJOTk4mMDDQpH2Zso9Tp06RlJTE888/n+uYz507h1arpUmTJoZlrq6uNGjQgFOnThmWPT4Wl5+fn2GW26zcv38/yz5B5cqVY/Xq1bRo0YKnn36a5cuXZ9qPJa1ASouhYcOGBAQE8O233xq9l5O0ifBM7YiaW3kuRrZv357ta4Bu3brRrVu3vB5KCCFU6auvYNky/fMXXni0fMeOnIsRNTh16hTly5cnPj4egN9++40yZcoYrWPqJQJT9pH2BWsLj/eZ0Gg0Od7+XKxYMe7evZvpe9HR0bz55pt07tyZ/fv38/777zN37twc4/D19eXZZ58lMjKSYsWK4ezsbFJ/zjt37gBQvHjxHI+RFzI3jRBC2NmuXfp/XVwg/aX/TP5vpzp//PEHx44d45VXXqFq1aq4u7tz6dIlKlSoYPTw9/c3aX+m7KNixYoUKFCAbdu2ZbkfNzc3UlMzjvWS5plnnsHNzY3du3cblmm1Wvbv30/VqlVNzD5zgYGBnDx5MsPyW7du8fzzz1OlShV+/PFHtm3bxjfffMMHH3yQ4z7j4+M5d+4cfn5+uLm5UbduXaP8dTodf/zxR4b+nMePH6ds2bIUK1YsTznlxOKXaYQQQpju1i142G+SunWhcWN4+WVo2BCy6VaQLyUlJREVFWW4tXfz5s2EhobywgsvEBwcjLOzMx988AHvv/8+Op2Opk2bEhMTw+7du/H29jb078iOl5dXjvvw8PDgww8/ZOTIkbi5udGkSRNu3rzJiRMn6N+/PwBPPfUUf//9NxcuXMDT0zPDbccFCxZk8ODBjBgxgiJFilCuXDk+++wzEhMTDfvIrXbt2jF69Gju3r1rGOVcp9PRoUMHAgIC+Oabb3BxcaFq1aqEh4fTqlUrypQpw/tpt2IBH3zwAZ07dyYgIIBr167x8ccf4+zsTI8ePQAICQmhd+/e1KtXj3r16vH555+TkJBA3759jWL566+/aNu2bZ7yMYUUI0IIYUd79jx63rQpODvDDz/YLx5r2rx5M35+foZBz2rVqsWcOXPo3bu3oa/h5MmTKV68OKGhofz333/4+vpSp04dPvroI5OPY8o+xo0bh4uLC+PHj+fatWv4+fkxaNAgw/sffPABvXv3pmrVqty/f5+jR49mOM7UqVPR6XS88cYbxMXFUa9ePbZs2ZLjNCk5qVGjBnXq1OHbb7/lrbfeAvR3u3z66ac0a9YMNzc3w7q1atXi999/z3AZ5cqVK/To0YPbt29TvHhxmjZtyt69ew3rvf7669y8eZPx48cTFRVlmMQ2fafWBw8esGHDBjZv3pynfEyi5AMxMTEKoMTExFh0v8nJycqGDRuU5ORki+7XUag9P0VRf45qz09R1J9jTvmNGKEo+vtmFGX9+uz3df/+feXkyZPK/fv3LR9oHqSmpip3795VUlNT7R2KVdgjv19//VWpUqWKTY6ZVX4LFixQ2rRpk+P22f1emvr9LS0jQghhR2n9RUB/iUYIgE6dOnH27FmuXr1qcn8ZS3N1dTWpc6wlSDEihBB2cv8+pN2Z++yz8HAUbgC0Wjh8GPbu1Q9+ZoNRyIWDST9YqD0MGDDAZseSYkQIIezkwAF90QH6/iLpdesGP/2kf96+vb5YEUKt5NZeIYSwk/SXaB4vRho1evT8SbjFVzzZpBgRQgg7qVkTuneHMmUyFiMPp1AB9IOfCaFmcplGCCHspFMn/UNRMr5Xrx54ekJi4qN5aoRQK2kZEUIIO9NoMnZQdXWFtGlPrl6F//6zfVxC2IoUI0II4aAaNnz0PG2UViHUSIoRIYSwg/v3IYf50kg/i/xjc5qJHFy4cAGNRsORI0fsHYowgRQjQghhB336gLu7vuC4ejXzddJPoBoVZZOwrKZPnz507drV3mEIByXFiBBC2EF0NKSkwPXr4Oub+TrppgmRlhGhalKMCCGEHaQVFwUL6h+ZKVkSChWCZ54BHx/bxWZrO3bsoEGDBri7u+Pn58eoUaNISUkxvL9582aaNm2Kr68vRYsW5YUXXuDcuXNG+9i3bx+BgYF4eHhQr149Dh8+bOs0RB5IMSKEEHZw44b+3/StH497+mmIi4PISJg0yTZx2drVq1fp2LEj9evX5+jRoyxcuJDly5fzySefGNZJSEggJCSEAwcOsG3bNpycnHjppZfQPex0Ex8fzwsvvEDVqlU5ePAgEyZM4IMPPrBXSiIXZJwRIYSwMa0W7tzRP8+uGDF1PpoZM/SPnNSpAz//bLzsxRfh0KGctw0J0T8sbcGCBfj7+zNv3jw0Gg2VK1fm2rVrfPjhh4wfPx4nJydeeeUVo21WrFhB8eLFOXnyJNWrV2ft2rXodDqWL1+Oh4cH1apV48qVKwwePNjyAQurkGJECCFsLK1VBLIvRkwVG5t1J9j0Mpv89eZN07aNjTU/LlOcOnWKoKAgNOkqryZNmhAfH8+VK1coV64cZ8+eZfz48fz999/cunXL0CJy6dIlqlevzqlTp6hZsyYeHh6GfQQFBVknYGEVUowIIYSNpe+Mmn6m3tzy9tYPKZ+T4sUzX2bKtt7e5sdlKZ07dyYgIIClS5dSunRpdDod1atXJzk52X5BCYuSYkQIIWzMnJaRxYthyxb9rb3ffJP5Onm5hPL4ZRtbq1KlCj/88AOKohhaR3bv3o2Xlxdly5bl9u3bnD59mqVLl9KsWTMAdqWfYfDhPlavXs2DBw8MrSN79+61bSIiT6QDqxBC2Fj6lpGcipHDh2H9eoiIMC5i8qOYmBiOHDli9HjzzTe5fPky77zzDv/++y8//fQTH3/8MSEhITg5OVG4cGGKFi3KkiVLiIyM5I8//iDkscqrZ8+eaDQaBg4cyMmTJ9m4cSPTp0+3U5YiN6RlRAghbMycYiT9+7duQbly1onJFrZv305gYKDRsv79+7Nx40ZGjBhBrVq1KFKkCP3792fs2LEAODk5sW7dOt59912qV69OpUqVmDNnDi1btjTso1ChQvzyyy8MGjSIwMBAqlatyrRp0zJ0fBWOS4oRIYSwseBg/bwz0dGQUz/L9MXI7dv5txgJCwsjLCwsy/f37duX5XutW7fm5MmTRsuUx6YxbtSoUYah3x9fRzguKUaEEMLGSpUyHuo9O4+3jAihRtJnRAghHNjjLSNCqJEUI0II4cDSt6BIMSLUSi7TCCGEjX31lX6umXLloGbN7NeVyzTiSSDFiBBC2FBqqr4Dq04HdevCgQPZr1+oEBQoAPfv60dLBemYKRyLJX4f5TKNEELY0O3b+kIETBsKXqN5tN7Zs64AJCYmWik6IcyX9vvo6uqa631Iy4gQQtiQOWOMpOndWz97r5+fMz4+vtx4OPqZp6en0Zwu9qLT6UhOTubBgwc4Oanv/7iSX+YURSExMZEbN27g6+uLs7NzrmMwqxhZuHAhCxcu5MKFCwBUq1aN8ePH06FDh0zXDwsLo2/fvkbL3N3defDgQe6iFUKIfC43k+RNmPDouaKUQqPBUJA4AkVRuH//PgUKFHCI4sjSJL/s+fr6UsrUe9WzYFYxUrZsWaZOnUrFihVRFIVVq1bRpUsXDh8+TLVq1TLdxtvbm9OnTxteq/FECiGEqfI6SZ5Go8HPz48SJUqg1WotF1geaLVadu7cSfPmzfPUVO+oJL+subq65qlFJI1ZxUjnzp2NXk+ZMoWFCxeyd+/eLIsRjUaT54pJCCHUIjeXaTLj7OxskS8BS3B2diYlJQUPDw9VfllLftaX6z4jqampfPfddyQkJBCUzXjG8fHxBAQEoNPpqFOnDp9++mmWhUuapKQkkpKSDK9jY2MBffVmyf8JpO3LUf53YWlqzw/Un6Pa8wP15/h4fteuOQH6IqJo0RS0WtPuRNDp9J1fXV3B19cakebek3YO1caa+Zm6T41i5j05x44dIygoiAcPHlCoUCHWrl1Lx44dM103IiKCs2fPUrNmTWJiYpg+fTo7d+7kxIkTlC1bNstjTJgwgYkTJ2ZYvnbtWjw9Pc0JVwghHMqcOYH88Yd+gpnZs/8gICAux2327SvJ1KkN0OmcCA4+wcsvR1o7TCEsIjExkZ49exITE4O3t3eW65ldjCQnJ3Pp0iViYmL4/vvvWbZsGTt27KBq1ao5bqvVaqlSpQo9evRg8uTJWa6XWcuIv78/t27dyjYZc2m1WsLDw2nTpo0qm97Unh+oP0e15wfqz/Hx/F580ZnNm/V3LFy5ojWp38iuXRpatdI3ZA8blsrnn+usGbLZnrRzqDbWzC82NpZixYrlWIyYfZnGzc2NChUqAFC3bl3279/P7NmzWbx4cY7burq6EhgYSGRk9lW9u7s77u7umW5vjV8Ea+3XUag9P1B/jmrPD9SfY1p+BQroR1+Nj4dSpVwxpdtHmTKPnt+86Yyrq2P0FXnck3IO1coa+Zm6vzzfMK3T6YxaMbKTmprKsWPH8PPzy+thhRAiX/rxR7h3DxISMKkQAeOOruk7wAqhFma1jIwePZoOHTpQrlw54uLiWLt2Ldu3b2fLli0ABAcHU6ZMGUJDQwGYNGkSjRo1okKFCty7d4/PP/+cixcvMmDAAMtnIoQQ+Ugmjb9Z8vbWr5+UJMWIUCezipEbN24QHBzM9evX8fHxoWbNmmzZsoU2bdoAcOnSJaPR2+7evcvAgQOJioqicOHC1K1blz179pjUv0QIIYRe2pDwly5JMSLUyaxiZPny5dm+v337dqPXM2fOZObMmWYHJYQQwlhaMXLrFqSkgItM5iFURH2D7AshhIM6fBhefBEGDoRNm8zbNm3sSEV5NHuvEGohtbUQQthIZCT88ov+ecWKkMW0Xpl6vBOr3Acg1ERaRoQQwkbyMhS83FEj1ExaRoQQwkbyUoz06QNt2+q3CwiwaFhC2J0UI0IIYSN5KUYqVNA/hFAjuUwjhBA2cuPGo+d5mbFXCLWRYkQIIWwkfctI8eL2i0MIRyOXaYQQwkbSipEiRcDcKUB0Oti4Ub+PAgWgZ0/LxyeEvUgxIoQQNpJWjOTmEo1GA6+8AsnJUKOGFCNCXeQyjRBC2EBCAiQm6p/nthhJG/hMbu0VaiMtI0IIYQOpqTBihL6QqFw5d/tIPyR8aqrps/4K4eikGBFCCBvw9obPPsvbPtJaVHQ6fUEid+QItZDLNEIIkU+kLz6iouwXhxCWJsWIEELkEzIkvFArKUaEEMIGEhP1/TzyQooRoVZSjAghhA2MHu2EqyuUKAHHjuVuH1KMCLWSYkQIIWwgOlqDosDNm/rOrLmRdmuvfn+WiUsIRyB30wghhA2kn5emRInc7aNUKfD11beQ+PpaIiohHIMUI0IIYQPR0RoAvLz0w7nnRqVKcPeuBYMSwkHIZRohhLCBmzf1/8rYIEJkJMWIEEJYmVbrxL17+paR3F6iEULNpBgRQggri4t7NEVvsWJ2DEQIByV9RoQQwsoSEtwMz4sUydu+Jk+Gffv0fUd27JD5aYQ6SDEihBBWlr5lpHDhvO1rzx7YvFn/PCYm78WNEI5ALtMIIYSVxcdbrhhJX3zcuZO3fQnhKKRlRAghrKxq1Tv88UcKsbEuPPts3vaVvpiR23yFWkgxIoQQVlaokJamTRVcXXNeNydSjAg1kss0QgiRj6QvRuQyjVALKUaEECIfSd9nRFpGhFrIZRohhLCyEyeKkJqqoXhxaNAA3N1zvy+5TCPUSIoRIYSwsh9/rMjBg/o/tzduQPHiud+XXKYRamTWZZqFCxdSs2ZNvL298fb2JigoiE2bNmW7zXfffUflypXx8PCgRo0abNy4MU8BCyFEfhMf/2jQs7zOtistI0KNzCpGypYty9SpUzl48CAHDhygVatWdOnShRMnTmS6/p49e+jRowf9+/fn8OHDdO3ala5du3L8+HGLBC+EEPlB2jgjXl7k+Y4aPz8YMABGjoROnSwQnBAOwKzLNJ07dzZ6PWXKFBYuXMjevXupVq1ahvVnz55N+/btGTFiBACTJ08mPDycefPmsWjRojyELYQQ+UdaMZLXAc9AP7fN0qV5348QjiTXfUZSU1P57rvvSEhIICgoKNN1IiIiCAkJMVrWrl07NmzYkO2+k5KSSEpKMryOjY0FQKvVotVqcxtyBmn7suQ+HYna8wP156j2/ED9OSYna0lI8ADA11dBq02xc0SWp/ZzKPnlfd85MbsYOXbsGEFBQTx48IBChQqxfv16qlatmum6UVFRlCxZ0mhZyZIliYqKyvYYoaGhTJw4McPyrVu34unpaW7IOQoPD7f4Ph2J2vMD9eeo9vxAvTk+eOBMSsoLAOh0t9i4cY+dI7IetZ7DNJKf+RITE01az+xipFKlShw5coSYmBi+//57evfuzY4dO7IsSHJj9OjRRi0qsbGx+Pv707ZtW7y9vS12HK1WS3h4OG3atMHVEkMjOhi15wfqz1Ht+YH6czx//lFLSIUKRenYsaNF9puSAvfu6S/b2Jvaz6Hkl3tpVzZyYnYx4ubmRoUKFQCoW7cu+/fvZ/bs2SxevDjDuqVKlSI6OtpoWXR0NKVKlcr2GO7u7rhnciO+q6urVX4RrLVfR6H2/ED9Oao9P1BvjnFxj54XLeqEq2vex5ps2xbS/hOblARubtmvbytqPYdpJL/c7dMUef5U6HQ6o/4d6QUFBbFt2zajZeHh4Vn2MRFCCLW5d09jeG6JDqxgPGia3N4r1MCslpHRo0fToUMHypUrR1xcHGvXrmX79u1s2bIFgODgYMqUKUNoaCgAw4YNo0WLFnzxxRd06tSJdevWceDAAZYsWWL5TIQQwgElJUGhQskkJLhSuLAm5w1M8PhYI491zRMi3zGrGLlx4wbBwcFcv34dHx8fatasyZYtW2jTpg0Aly5dwsnpUWNL48aNWbt2LWPHjuWjjz6iYsWKbNiwgerVq1s2CyGEcFBt2iisWbOJ9u07otFYpglc5qcRamNWMbJ8+fJs39++fXuGZd26daNbt25mBSWEEGrj5JT3Ac/SyCisQm1k1l4hhMhnZH4aoTZSjAghRD4jl2mE2kgxIoQQVjRrlhNz59Zm5Egn7t2zzD7lMo1QGylGhBDCin7/XcO2bQHMmuWMTmeZfUoxItRGihEhhLCitGJBo1Hw8bHMPtNfppE+I0INcj1RnhBCiJzdvasfW8THB5ydLbPPcuXg11/1RUnZspbZpxD2JMWIEEJYUVrLiKVGXwXw9IROnSy3PyHsTS7TCCGElSgKhk6rhQsrdo1FCEcmxYgQQlhJXBykpuov01iyZUQItZHLNEIIYSXpO5f6+lp23/v3w+XL+paXfv0su28hbE2KESGEsJL0t91a+jLNiBGwY4f+eY8eUKCARXcvhE3JZRohhLAS42LEsvuWsUaEmkgxIoQQVlKsGPTuraNhw+vUqGHZlhEZEl6oiVymEUIIK6lZE5YuTWXjxn107NjRovuWyfKEmkjLiBBC5ENymUaoiRQjQgiRD8mQ8EJNpBgRQggrUaw4zpm0jAg1kWJECCGspGdPKFnShbfeas3165bdtxQjQk2kA6sQQljJrVtpE+UVpGBBrUX3LZdphJpIMSKEEFaS1mLh5KTDy8uy+y5cGFxc9P+6u1t230LYmhQjQghhJWnFSMGCWjQay14Vf+YZSE4GjcaiuxXCLqTPiBBCWEna5RMvL8teogF9ESKFiFALKUaEEMIKdDqIidE/t3R/ESHURooRIYSwgpiYR7f2FiqUbN9ghHBw0mdECCGsIP3ttoUKWadl5Isv4OhRuHcPfvpJLtuI/EuKESGEsIL0t9taqxjZtAm2bdM/j4/H4nfsCGErcplGCCGswLhlxDqXaWSyPKEW0jIihBBWULu2/tLJzZsp3Lt3HXja4sdIP/DZ3bsQEGDxQwhhE1KMCCGEFRQvDi++CFqtwsaNMVY5hgwJL9RCLtMIIUQ+JZdphFpIMSKEEPmUtIwItTCrGAkNDaV+/fp4eXlRokQJunbtyunTp7PdJiwsDI1GY/Tw8PDIU9BCCOHojhyBnTvh+HHQaq1zz+3jfUaEyK/MKkZ27NjBkCFD2Lt3L+Hh4Wi1Wtq2bUtCQkK223l7e3P9+nXD4+LFi3kKWgghHN3kydCiBdSp48q9e9aZyU4u0wi1MKsD6+bNm41eh4WFUaJECQ4ePEjz5s2z3E6j0VCqVKncRSiEEPmQLQY9k8s0Qi3ydDdNzMOJF4qkbyvMRHx8PAEBAeh0OurUqcOnn35KtWrVslw/KSmJpKQkw+vY2FgAtFotWq3lPtRp+7LkPh2J2vMD9eeo9vxAvTneueMCaHB1VfDwSLVKfsWKwcsvO1O4MAQF6dBqFYsfwxRqPYdpJL+87zsnGkVRcvXbq9PpePHFF7l37x67du3Kcr2IiAjOnj1LzZo1iYmJYfr06ezcuZMTJ05QtmzZTLeZMGECEydOzLB87dq1eHp65iZcIYSwqYED23Dzpic+Pg9YtWqLvcMRwi4SExPp2bMnMTExeHt7Z7lerouRwYMHs2nTJnbt2pVlUZEZrVZLlSpV6NGjB5MnT850ncxaRvz9/bl161a2yZhLq9USHh5OmzZtcHV1tdh+HYXa8wP156j2/EC9ORYt6kJcnIZnn9Xx2We/qC6/9NR6DtNIfrkXGxtLsWLFcixGcnWZZujQofz666/s3LnTrEIEwNXVlcDAQCIjI7Ncx93dHXf3jB2+XF1drfKLYK39Ogq15wfqz1Ht+YG6ckxJgbg4/fO0q9hqyi8ras9R8svdPk1h1t00iqIwdOhQ1q9fzx9//EH58uXNDiw1NZVjx47h5+dn9rZCCJEf3Lv36Hn6TqbWlJJim+MIYQ1mFSNDhgxhzZo1rF27Fi8vL6KiooiKiuL+/fuGdYKDgxk9erTh9aRJk9i6dSv//fcfhw4d4n//+x8XL15kwIABlstCCCEcSPrbbK1djLz4on62Xnd3SE217rGEsBazLtMsXLgQgJYtWxotX7lyJX369AHg0qVLODk9qnHu3r3LwIEDiYqKonDhwtStW5c9e/ZQtWrVvEUuhBAOKv1ttoULW/cOl+RkiI/XP4+JMR4ITYj8wqxixJS+rtu3bzd6PXPmTGbOnGlWUEIIkZ/FxYFGA4oCvr7WPdbjY41IMSLyI5m1VwghLKx1a30fjthYSE3Vkc3oB3kmQ8ILNZBiRAghrMDJSd8qYu1xsmRIeKEGMmuvEELkYzIkvFADKUaEECIfS3+ZRlpGRH4ll2mEEMLCVqyA48f1hcKbb1r3WNIyItRAihEhhLCwn3+Gn37SPw8Otu6xpBgRaiCXaYQQwsJsOeiZXKYRaiAtI0IIYWFpLRQFCoCHh3WP9dRTsGaNvuh55hnrHksIa5FiRAghLCythcIW89J4eUGvXtY/jhDWJJdphBDCwtJaRmw1SZ4Q+Z0UI0IIYUEPHkDa3KEyNLsQppHLNEIIYUHGk+TZ5pinTsH16/qJ8l56yTbHFMKSpBgRQggLSn9Hi61aRgYMgD179M+Tk8HV1TbHFcJS5DKNEEJYkD1aRmSsEZHfScuIEEJYUKFC0LWrviioUsU2x3y8GClRwjbHFcJSpBgRQggLql0b1q9/9Nras/aCDHwm8j+5TCOEEPmcXKYR+Z0UI0IIkc+lbxmRYkTkR1KMCCFEPpe+ZUQu04j8SIoRIYSwoEGD4OmnoV49uHrVNseUyzQiv5MOrEIIYUGXLsH58/qHu7ttjikdWEV+Jy0jQghhQelbJnx9bXPMtJYRV1dISbHNMYWwJGkZEUIIC0orRry9wcXFNrf2PvssxMVBwYKg0Vj/eEJYmhQjQghhQWmXSWw5Y6+zs36wNSHyK7lMI4QQFqIoj1pGZMZeIUwnxYgQQlhIfPyjPhu2bBkRIr+TyzRCCGEh6Tuv2rplZNEiOHlS33dk5UrbHluIvJJiRAghLMQeM/am+fpr2LlT/3zBAihQwLbHFyIv5DKNEEJYSPoxPmzdMiIDn4n8TFpGhBDCQipXhrAwfTFQr55tj/34wGelS9v2+ELkhRQjQghhIX5+0Lu3fY4tLSMiPzPrMk1oaCj169fHy8uLEiVK0LVrV06fPp3jdt999x2VK1fGw8ODGjVqsHHjxlwHLIQQIiOZuVfkZ2YVIzt27GDIkCHs3buX8PBwtFotbdu2JSEhIctt9uzZQ48ePejfvz+HDx+ma9eudO3alePHj+c5eCGEEHoyc6/Iz8y6TLN582aj12FhYZQoUYKDBw/SvHnzTLeZPXs27du3Z8SIEQBMnjyZ8PBw5s2bx6JFi3IZthBCOJ6zZ+H+fX0rhZ+ffmRUW5HLNCI/y1OfkZiYGACKZNNtPCIigpCQEKNl7dq1Y8OGDVluk5SURFJSkuF1bGwsAFqtFq0FJ3pI25cl9+lI1J4fqD9HtecH6spx9GhnfvhB3+B85oyWp56yXX7e3hrS/qTfvJmKVquz6vHSU9M5zIzkl/d95yTXxYhOp+O9996jSZMmVK9ePcv1oqKiKFmypNGykiVLEhUVleU2oaGhTJw4McPyrVu34unpmduQsxQeHm7xfToStecH6s9R7fmBOnKMjAwCSgCwf/9WTp58NIWutfM7c8YXaAHA0aMX2bjxmFWPlxk1nMPsSH7mS0xMNGm9XBcjQ4YM4fjx4+zatSu3u8jS6NGjjVpTYmNj8ff3p23btnh7e1vsOFqtlvDwcNq0aYOrq6vF9uso1J4fqD9HtecH6spx4kT9n1RnZ4VXX22LRmO7/KpXh82bdfj6Qrt25ejY0d9qx3qcms5hZiS/3Eu7spGTXBUjQ4cO5ddff2Xnzp2ULVs223VLlSpFdHS00bLo6GhKlSqV5Tbu7u64u7tnWO7q6mqVXwRr7ddRqD0/UH+Oas8P1JHjvXv6f319Nbi5Gedi7fyeeQa2bUt7ZZ/xLNVwDrMj+eVun6Yw6zdWURSGDh3K+vXr+eOPPyhfvnyO2wQFBbHt0ScE0DcFBQUFmXNoIYRweDJjrxC5Y1bLyJAhQ1i7di0//fQTXl5ehn4fPj4+FHg4EUJwcDBlypQhNDQUgGHDhtGiRQu++OILOnXqxLp16zhw4ABLliyxcCpCCGE/qamPWkZkxl4hzGNWy8jChQuJiYmhZcuW+Pn5GR7ffPONYZ1Lly5x/fp1w+vGjRuzdu1alixZQq1atfj+++/ZsGFDtp1ehRAiv3l4cyEgLSNCmMuslhFFUXJcZ/v27RmWdevWjW7duplzKCGEyFfsOWNvmh49YPduiI+H27dBo7FPHEKYS+amEUIIC7DnjL1pbt6Ey5f1z+PjwcvLPnEIYS77dLkWQgiVcYSWERkSXuRX0jIihBAW0Lq1vgPrnTtQsKB9Ynh8sryAAPvEIYS5pBgRQggLcHICHx/9w16kZUTkV3KZRgghHISiKNnOgp6Tx1tGhMgvpBgRQgg7u337NtOnT6dixYp4eXkxcOBA4uPjzd6PzNwr8iu5TCOEEBawbh2cPatvnfjf/0y7XHPgwAHmzZvHunXrjGYqX7ZsGdu3b+err76iQYMGJscgl2lEfiUtI0IIYQHffAPjx8PQofrbanMyZcoU6tevz6pVq4wKkbR5uSIjI2ncuDFTpkwhNTXVpBjkMo3Ir6QYEUIICzB1nBFFURgzZgxjx441LPP19eX999/n9OnTHD9+3NAakpqaytixY2nVqhV3Tagu5DKNyK+kGBFCCAtI+/J3d4eHU3VloCgKISEhfPrpp4ZlEydO5OrVq8yYMYNnn32WChUqsGvXLsaNG4eTk/5P9M6dO2ndujW3b9/ONoannoK5c2HNGnjrLUtkJYRtSDEihBAWkNOMvTqdjsGDBzNr1izDsnnz5jF+/Hg8PT2N1nV1dWXSpEns2LGDEiVKAHDo0CFatWrFjRs3soyhcGH9ZaJevSAwME/pCGFTUowIIYQFpF2myWz0VUVRGDx4MIsXLwZAo9GwfPlyhgwZku0+mzZtyvbt2/Hz8wPgn3/+4bnnnjOajFQINZBiRAgh8ig5GRIT9c8zaxnZtGkTK1euBMDZ2ZmvvvqKfv36mbTvKlWqsGPHDsqWLQvAyZMnadmyJVevXrVI7EI4AilGhBAij7Kbl2b//v2sWLHC8HrNmjX06NHDrP1XrFiRHTt2EPBwfPczZ87QrVs3dDpdhnWvXIFDhyCTCdSFcFhSjAghRB5ldSfN7du36dGjBykpKQCEhITQvXv3XB3j6aefNipIIiIiWLRoUYb1Xn0V6taF554DE+8IFsLupBgRQog8yqxlRKfT8cYbb3Dp0iUAgoKCmDp1ap6OExAQQFhYmOH1qFGjMlyuSd8yExOTp8MJYTNSjAghRB65u0PLllCzJpQvr18WGhrKpk2bAPDx8WHt2rW4urrm+VgtW7akf//+AMTFxfHOO+8YvS9jjYj8SIoRIYTIo7p14c8/4ehRePdd/bgg48ePB/R3zoSEhFCmTBmLHe+zzz4z3PK7fv161q9fb3gv/WUiGRJe5BdSjAghhAWlpKQwZMgQQ+fS8ePHU6tWLYseo0iRIsyePdvweujQocQ8vCYjLSMiP5JiRAghLGjZsmUcP34cgPr16zN69GirHOf111+nY8eOAFy7do2PPvoIkGJE5E9SjAghhIXcu3ePcePGGV7PmjXLMKS7pWk0GhYsWGAYvXXhwoVERETIZRqRL0kxIoQQefTBB/p+I7Vq3eDWLRcAunfvTuPGja163ICAAD755BNAP8rrm2++SaFCWsP70jIi8gspRoQQIo9OndIPNHbp0rNACh4eHnm+jddU77zzDnXq1AHg+PHjbNmyzvCeFCMiv5BiRAgh8sj4S/8eH3zwgWFwMmtzcXFh6dKlhstBq1bNfBTJPZuEIESeSTEihBB5dOVKwsNnsfj5FefDDz+06fHr1KnDe++9B4BWe4wmTboTF6ewZIlNwxAi16QYEUKIPEhNTeX69QcPX90lNDSUQoUK2TyOiRMnUq5cOSCF3bu/Yf36NWg0Ng9DiFyRYkQIIfJg+fIVpKR4AVCgwAPeeOMNu8RRqFAhFixYYHj9/vvvc+vWLbvEIoS5pBgRQohcio2NZcyYTwE3ACpXLmm1W3lN0alTJ7p16wboJ+mz1hgnQliaFCNCCJFLU6dO5dYtxfD6mWd87RfMQ7Nnz8bTMxiYwPLldYiM/M/eIQmRIylGhBAiFy5cuMCMGTOA4oZlD6eLsSs/Pz9Klx4HfIyiDGbixDn2DkmIHEkxIoQQuTB69GiSkpJwtGIEoE6dcobnX3/9OxcvXrRjNELkzOxiZOfOnXTu3JnSpUuj0WjYsGFDtutv374djUaT4REVFZXbmIUQwq4iIiJYt04/uFjhwtFMm3afsWOhZUv7xpWmbFk3w/PU1MI2G4BNiNwyuxhJSEigVq1azJ8/36ztTp8+zfXr1w2PEo7yXwghhDCDoiiEhIQYXk+ZMoCRIwsweTK0aGHHwNIpXjz9qxIsX76cy5cv2yscIXLkYu4GHTp0oEOHDmYfqESJEvj6+pq9nRBCOJJ169axd+9eAKpUqcLAgQPtHFFGxsVIcbRaLdOmTWPevHn2CkmIbJldjORW7dq1SUpKonr16kyYMIEmTZpkuW5SUtLDa7F6sbGxAGi1WrRabVabmS1tX5bcpyNRe36g/hzVnh/krxzj4+MZMWKE4fW0adNQFCXb2O2RX5EiGtL+vLu6lkGrhaVLl/LBBx9QpkwZix8vP53D3JD88r7vnGgURVFyXi2LjTUa1q9fT9euXbNc5/Tp02zfvp169eqRlJTEsmXLWL16NX///bdhcqfHTZgwgYkTJ2ZYvnbtWsN02UIIYWurV6/mhx9+AKBevXqMHTuW69cL4uKiw9s7CXd3nZ0j1Dt9ujAfftgcgKef3sh//3UC4IUXXmDAgAH2DE08YRITE+nZsycxMTF4e3tnuZ7Vi5HMtGjRgnLlyrF69epM38+sZcTf359bt25lm4y5tFot4eHhtGnTBldXV4vt11GoPT9Qf45qzw/yT45nzpwhMDAQrVaLm5sbR44coUKFCgQFOXPwoBNOTgoJCSk4OxtvZ4/8zp2DKlX0x3rxxUTCw4tx//59PDw8OHPmDKVKlbLo8fLLOcwtyS/3YmNjKVasWI7FiM0u06TXoEEDdu3aleX77u7uuLu7Z1ju6upqlV8Ea+3XUag9P1B/jmrPDxw7R0VRGD58uKHJecSIEVSpUgWAmzf16xQtqsHDI+v4bZlf6dKPnsfHezJ48GBmzJjBgwcPmD9/vtXurnHkc2gJkl/u9mkKu4wzcuTIEfz8/OxxaCGEMNtPP/3Eli1bAChXrhwfffQRAIryqBhxpBsEvbygfn1o0wYaNIDhw4fj5qa/3XfhwoXExMTYOUIhjJndMhIfH09kZKTh9fnz5zly5AhFihShXLlyjB49mqtXr/Lll18CMGvWLMqXL0+1atV48OABy5Yt448//mDr1q2Wy0IIIazk/v37vPfee4bXM2bMMPRdS0iA+/f1y43vYLEvjQb27Uu/pDS9e/dm6dKlxMbGsnDhQkaNGmWv8ITIwOyWkQMHDhAYGEhgYCAAISEhBAYGMn78eACuX7/OpUuXDOsnJyczfPhwatSoQYsWLTh69Ci///47zz//vIVSEEII65k2bZphBNPWrVvz8ssvG95LaxUBx2oZyczIkSMNk/jNnDmT+2lVlBAOwOyWkZYtW5Jdn9ewsDCj1yNHjmTkyJFmByaEEPZ24cIFpk2bBoCLiwtz5sxBo9EY3r9x49G6jtQykpkKFSrw6quv8u2333Ljxg3CwsIYPHiwvcMSApC5aYQQIkvDhw/nwYMHAAwbNszQaTVNfmoZAYwuzXz++eekpKTYMRohHpFiRAghMrFt2zZ+/PFHAEqWLGm4FJ1e+mLE0VpGFi2CWrX0d9bs2aNfFhgYSLt27QB9f79vv/3WjhEK8YgUI0II8ZiUlBSGDRtmeD116tRMx0hw5Ms0d+/CP//A9euQfl7S0aNHG55PnTo128vuQtiKFCNCCPGYhQsXcuLECUA/LlJwcHCm6znyZZr0xVH6OJs3b06jRo0AOHbsGBs3brRxZEJkJMWIEEKkc/PmTaNLMnPnzjXchfK4Tz6BixfhwAF4eIOhw0hfHKUvRjQajVHfkUmTJknriLA7KUaEECKdsWPHcu/ePQD69OlDgwYNslzXwwPKlYO6dfUDjTmSrFpGADp37kyNGjUA2Ldvn2G+HSHsRYoRIYR46PDhwyxduhQALy8vQkND7RxR7mVXjDg5ORkNCT969GiSk5NtFJkQGUkxIoQQQGpqKm+99ZbhksX48eMtPqGcLaW/TJO+o22aDh068NxzzwEQGRnJkiVLbBSZEBlJMSKEEMD8+fPZv38/AFWqVOHdd9/Ndn1FgQkTYO5ceDhtjUPx8oKH09FkaBkBfd+Rzz77zPB64sSJxMbG2ig6IYxJMSKEeOJdunTJMPkdwNKlSw0Ty2UlPh4mToR334V03+kOQ6N5dKkms2IEoF69evTo0QOAW7duGUabFcLWpBgRQjzRFEVhyJAhJCQkADBo0CCaNGmS43aOPMZImrS4bt3St+RkZsqUKYZp3mfMmMGVK1dsFJ0Qj0gxIoR4on3//ff8+uuvAPj5+ZncadWRR19NExIC8+bBmjWg02W+Tvny5Rk6dCgADx484OOPP7ZhhELoSTEihHhi3b17l3feecfweu7cufj6+pq0rSMPeJbmjTdgyBB47TVwds56vTFjxuDj4wPAypUrOXr0qI0iFEJPihEhxBNr1KhRREdHA/Diiy/y8ssvm7xtfrhMY6qiRYsyZswY4NFlKxkITdiSFCNCiCfSH3/8YbidtVChQsybNw+NRmPy9vnhMo053n33XSpWrAjA7t27Wb16tZ0jEk8SKUaEEE+c+Ph4BgwYYHgdGhqKv7+/WfvID5dpkpPhwgXYvx8uXcp+XXd3d+bOnWt4PWLECMNItEJYmxQjQognzkcffcT58+cBaNasGW+//bbZ+8gPl2k2bYLy5aFBAzCloaNdu3a88sorANy4ccNojh4hrEmKESHEE+Wvv/4ytAB4eHiwfPnyLCfCy05+aBlJXyRlNgprZmbOnImnpyegHwju8OHDVohMCGNSjAghnhiJiYn069fP8HrKlCmGfhLmqlBBP1NvuXJg4g04Npfd/DRZ8ff3N7SI6HQ6hgwZgi6r+4KFsBApRoQQT4xx48YRGRkJQFBQEMOGDcv1vubNg0OH4OJFyEXDik3kphgBeP/996lUqRIAERERhIWFWTYwIR7joB8hIYSwrN27dzNz5kxA31lzxYoVOGc3+IYK+PjAw8FVzSpG3NzcmDdvnuH1yJEjuWnODoQwkxQjQgjVu3DhAi+//LJh7IyJEydSuXJlO0dlfabMT5OV1q1b8/rrrwNw+/Zthg8fbuHohHhEihEhhKrFxMTQqVMnbjzswdmiRYsn6os1fTFi7jhms2bNMoxIu3r1asLDwy0bnBAPSTEihFAtrVZLt27dOHnyJADPPvssP/74Iy4uLnna76FDULMmtG4Nq1ZZIlLrSStGtFqIiTFv21KlSjF9+nTD60GDBpGYmGjB6ITQk2JECKFKiqIwdOhQw//mixYtym+//UaRIkXyvO8rV+DYMdi2DS5fzvPurCq3nVjT9OvXjxYtWgDw33//MXHiRAtFJsQjUowIIVTpiy++MAz37ubmxoYNG6hQoYJF9p0fBjxLkz6+27fN316j0bB48WLc3NwA/c/1yJEjlglOiIekGBFCqM769esZOXKk4fWKFSto2rSpxfafHwY8S/PRR/qWnAcPoFGj3O2jUqVKjB07FoDU1FQGDhxIamqqBaMUTzopRoQQqnLw4EF69epluHNmwoQJ9OrVy6LHyE+T5JUsCWXKgLt73vbz4YcfUrVqVQAOHDhg1JdEiLySYkQIoRpXrlyhc+fO3L9/H4BevXpZZX6V9JdpHL1lxFLc3NxYtmyZYej8cePGcfToUTtHJdRCihEhhCrEx8fTuXNnrl+/DkCTJk1YtmwZGo3G4sfKTy0jlhQUFGS4/KXVannjjTdISkqyc1RCDaQYEULke6mpqfTs2dPQsfLpp59m/fr1eHh4WOV4aS0jLi6OOy9Nmvh4mD0bxoyB5cvzvr8JEyZQs2ZNAI4dO8a4cePyvlPxxDO7GNm5cyedO3emdOnSaDQaNmzYkOM227dvp06dOri7u1OhQgWZ50AIYTE6nY63336bX375BQAfHx9+/fVXiluxySKtZaR4cf0op44sJQXeew8+/RS+/Tbv+3N3d2fNmjWGu2umT5/Orl278r5j8UQzuxhJSEigVq1azJ8/36T1z58/T6dOnXjuuec4cuQI7733HgMGDGDLli1mByuEEOnpdDoGDx5suIXX2dmZ77//nipVqljtmIryqBjJD/1F0s9Pk76vS17UqFGDTz75BNCP59K/f39DPx0hcsPsYQg7dOhAhw4dTF5/0aJFlC9fni+++AKAKlWqsGvXLmbOnEm7du3MPbwQQgAZCxEnJydWr15N69atrXxcmDVL/8VeuLBVD2URGg0UKwbXr+du0LOshISE8Msvv/DXX39x/vx5li1bxiuvvGK5A4gnSt7GRDZBREREhj8O7dq147333stym6SkJKNOUbGxsYC+w5RWq7VYbGn7suQ+HYna8wP156j2/CB3Oep0OoYOHcqyZcsAfSGyatUqXn31VZv8rAYMePQ8p8M5wjksXtyF69c13LypkJycYrFLS8uWLaNu3brEx8ezbds21q1bR/fu3S2zcwfiCOfQmqyZn6n7tHoxEhUVRcmSJY2WlSxZktjYWO7fv0+BAgUybBMaGprpkMNbt27F09PT4jGqffIntecH6s9R7fmB6TnqdDoWLVrE1q1bAX0h8v777+Pl5cXGjRutGWKe2PMcajRBQAmSkzX88MNWPD1TLLbvfv36MWfOHEA/d839+/cz/M1XC7V/Dq2Rn6lzGVm9GMmN0aNHExISYngdGxuLv78/bdu2xdvb22LH0Wq1hIeH06ZNG1zTLqqqiNrzA/XnqPb8wLwcFUXh3XffNSpEvvzyS1577TVbhJorjnAOv/7ambQhQWrXbouFRsUH9Jfur1+/znfffUdiYiJhYWFs27Ytz5MROhJHOIfWZM380q5s5MTqvy2lSpUiOjraaFl0dDTe3t6ZtoqAvre2eybDBbq6ulrlF8Fa+3UUas8P1J+j2vODnHNUFIVhw4axePFiQF+IrF27ltdff91WIQL6fhfx8fo7aQoWNP1uGnuew/QNFffuuWLpMBYsWMDOnTuJjo4mIiKC0NBQJk2aZNmDOAC1fw6tkZ+p+7P6OCNBQUFs27bNaFl4eDhBQUHWPrQQQiUUReGDDz5g7ty5wKPOqrYuRACWLoWnnwYvL/j5Z5sfPlfSFyNRUZbfv4+PDyEhITg7OwMwZcoUduzYYfkDCdUyuxiJj4/nyJEjhsGFzp8/z5EjR7h06RKgv8QSHBxsWH/QoEH8999/jBw5kn///ZcFCxbw7bff8v7771smAyGEqimKwujRo5kxYwagn0V25cqV9OzZ0y7xpG/oLVbMLiGYzd//0fOHf6otrlKlSkyYMAHQ9+v53//+x01L3r4jVM3sYuTAgQMEBgYSGBgI6G/vCgwMNMz/cP36dUNhAlC+fHl+++03wsPDqVWrFl988QXLli2T23qFEDlKSUlh2LBhTJs2zbBs2bJlRv/hsbXLlx89L1PGbmGY5ZlnoEYN6NQJypa13nE++OADnnvuOUA/T1DXrl158OCB9Q4oVMPsPiMtW7Y0zIaZmcxGV23ZsiWHDx8291BCiCdYbGws3bt3Z9OmTYZlixYtol+/fnaMCi5e1P/r5JR/ipHGjeGff6x/HGdnZ9asWUP9+vW5du0ae/bsoV+/fnz11VdWmSNIqIfMTSOEcDjnz5+ncePGhkLExcWFFStW8NZbb9k5skfFSJkyWLwjqBqULl2an3/+2TAMw9dff83HH39s56iEo5NiRAjhUPbs2UPDhg05ceIEAEWKFCE8PJy+ffvaOTJISIDbt/XPAwLsG4sjq1u3LmvXrjW0hkyePJkvv/zSzlEJRybFiBDCYaxatYrnnnvO0PHx2WefZe/evbRs2dK+gT2U1ioC+bcYURT9w9q6dOlimAYEYMCAAXKHjciSFCNCCLtLTU1l5MiR9OnTh+TkZABatWrF3r17qVixop2jeyQ/FyPTp+s7sfr4wJkztjnme++9x6BBgwD9wFovv/wyZ8+etc3BRb4ixYgQwq7u3bvHlClTmDVrlmHZ4MGD2bx5M4UdbCa6/FyM3L4Nx49DXJxxHtak0WiYO3eu4e7JO3fu8MILL3D37l3bBCDyDSlGhBB2c+rUKZo2bcqhQ4cAfUfVBQsWsGDBAocc6TL9bb1PPWW3MHIlffFkq2IE9Of0m2++oWrVqgCcOXPGZhMaivxDihEhhF2sWrWKevXqcebhNYOiRYsSHh7O4MGD7RxZ1j75RD/o2b590KiRvaMxT7lyj55ba+CzrPj4+PDrr79SvHhxAP744w/efvvtbIeJEE8WKUaEEDYVHx9PcHAwffr0MczoGRAQwJ49exymo2pWNBooUQLq1wcLztlpE/ZqGUlTvnx5NmzYYJh3bNmyZYZRdYWQYkQIYTNHjx6lXr16rF692rCsf//+fPbZZ5QvX96OkamfPVtG0jRu3JgVK1YYXo8YMYI1a9bYJxjhUKQYEUJYnU6nY/bs2TRs2JDTp08D4OXlxddff83ChQsznaVbWJaXF6T1B7ZHy0ianj17GqYPURSF4OBgGYNEmD8cvBBCmOP69ev07duXLVu2GJbVqVOHb775hgoVKuSbjozXr8O0afrLHY0bQ8OG9o7IfAEBcPcuXLkCqanwcJJdm5swYQK3b99m/vz5KIpCnz59UBSF3r172ycgYXfSMiKEsJqff/6ZmjVrGhUiw4cPZ8+ePVSoUMGOkZnv339h9mwICYHvvrN3NLmTdqkmJUVfXNlL2i2/Q4cOBfQtJH379s10bjPxZJBiRAhhcTExMQwYMIAuXbpw69YtAPz8/Ni6dSvTp0/Pl5dl8vMYI2ns3Yk1PY1Gw5w5c3jnnXcAfUHSr18/KUieUHKZRghhUVu2bGHAgAFcuXLFsKxLly4sW7aMYsWK2TGyvFFDMfLqq1Clij7+h8N+2JVGo2H27Nk4OTkxe/ZsQ0Hi6upKr1697B2esCEpRoQQFhEbG8vw4cNZtmyZYVmhQoX44osvGDhwYL6fQl4NxUjz5vqHI9FoNMycORPAUJAEBwfj5uZGt27d7BydsBUpRoQQebZhwwbeeecdo9aQ559/nuXLlxOQX7+5H6OGYsRRpRUkWq2WBQsWoNPp6NGjB66urnTt2tXe4QkbkD4jQohcu3TpEl26dOGll14yFCKFChVi4cKFhIeHq6YQAbhwQf+vtzf4+tozEnVK69Q6YMAAQD954muvvcZvv/1m58iELUgxIoQwW0pKCjNmzKBq1ar8/PPPhuVt27bl2LFjDBo0KN9flklPp3s0L01+r6+uXIFdu+CXX+wdSUZOTk4sXryY4OBg4NFMv+vXr7dzZMLapBgRQpjlr7/+ok6dOgwfPpyEhAQASpYsyddff83mzZt5Kr/NIGeC69chbTiU/F6MtGoFzZpBz57giFPDODk5sWLFCrp37w5AcnIyr776KosXL7ZzZMKapBgRQpgkKiqK4OBgmjdvzrFjxwB90/rgwYP5999/6d69u6paQ9JTU3+RtPjj4/UDoDkiZ2dnVq9ezRtvvAHoR/AdNGgQEyZMkMn1VEqKESFEtpKSkpg5cyaVKlUymlOmTp067NmzhwULFuCr8k4UBQrAK69A3bpQrZq9o8mb9MWUveaoMYWLiwthYWGMGDHCsGzixIkMGjSI1NRUO0YmrEHuphFCZCo1NZW1a9cyfvx4LqT13gQKFy7MlClTePPNN3G213jiNhYYCN9/b+8oLCP9hHkXL0Lt2nYLJUdOTk589tln+Pn5ERISAsCSJUu4ceMGX3/9NR4eHnaOUFiKtIwIIYwoisLGjRupU6cOwcHBRoXIgAEDOHPmDIMHD35iChG1yS8tI+m9//77fPXVV7i6ugL6W8nbt29PTEyMnSMTliLFiBDCYO/evbRs2ZJOnTrxzz//GJa3a9eOQ4cOsXTp0nw9iqrI2DKSX/Ts2ZNff/2VggULArBjxw5atGhBVFSUnSMTliDFiBCCU6dO8dJLLxEUFMTOnTsNy+vXr8+2bdvYvHkzgYGBdozQvtTUZ9KR5qcxV9u2bfnzzz8pWrQoAEePHqVJkyacO3fOzpGJvJJiRIgn2H///Ue/fv2oXr06GzZsMCyvWLEi3377LX///TetWrWyX4AOQFGgRAl49lno39/e0eRd2bKQdtNTfrlMk179+vXZvXs35R428fz33380adKEv//+286RibyQYkSIJ1BkZCR9+/bl2WefZeXKleh0OkA/s+6iRYs4ceIE3bp1U+2tuua4cwdu3YKzZx8NfJafubmBn5/+eX5rGUlTqVIl9uzZQ7WHtzZFR0fTvHlzVq1aZefIRG5JMSLEE+TMmTMEBwdTqVIlwsLCDLdI+vj4EBoaSmRkJG+99Zaho6BQ1xgjaQICwNkZPDweDeaW35QpU4adO3fSokULQD84Wp8+fRg+fDgpKSl2jk6YS4oRIZ4Ap0+f5o033qBKlSqsXr3a0BJSuHBhJk2axMWLFxk1ahSenp52jtTxqLEY+fVXePBAP99Ofq47ixQpQnh4OIMHDzYsmzFjBp06deKuo47oJjIlxYgQKnbq1Cn+97//UbVqVdasWWMoQooUKcInn3zChQsXGDduHD4+PnaO1HGpsRgpUgRcVDLKlKurKwsWLGDRokW4PExq69at1K9f3zBSsHB8UowIoTKKorB161Y6duxI1apV+eqrrwxFSNGiRfn00085f/48Y8aMwdvb287ROj41FiNq9NZbb7Ft2zbDrefnzp2jUaNGrFu3zs6RCVPkqhiZP38+Tz31FB4eHjRs2JB9+/ZluW5YWBgajcboIaPmCWF5iYmJLF68mGrVqtGuXTs2bdpkeK9o0aKEhoZy/vx5Ro8eLUWIGaQYyT+aN2/OwYMHqVu3LqD/TPTo0UP6keQDZhcj33zzDSEhIXz88cccOnSIWrVq0a5dO27cuJHlNt7e3ly/ft3wuJhfu3AL4YAuX77MqFGjKFu2LIMGDeLUqVOG9wICApg+fToXLlxg1KhReHl52THS/Cntz5WzM5QpY99YLOXuXRg9Wj9z76ef2jsayypXrhx//fUXffr0MSybMWMGbdq0kQHSHJjZxciMGTMYOHAgffv2pWrVqixatAhPT09WrFiR5TYajYZSpUoZHiVLlsxT0EI86RRFISIigu7du1O+fHmmTZtm1GGvWbNm/PDDD0RGRjJ8+HAKFSpkx2jzt7RipEwZ9fSzcHaGqVPh668hPNze0VhegQIFWLFiBQsWLDDcGbZ9+3YCAwPZsWOHnaMTmTHro5WcnMzBgwcZPXq0YZmTkxOtW7cmIiIiy+3i4+MJCAhAp9NRp04dPv30U8P94ZlJSkoiKSnJ8Do2NhYArVaL1oL3oaXty5L7dCRqzw/Un+Pj+cXFxbFu3ToWL15sNFw76Dvyvf766wwdOpQ6deoA+qLF0X82jnwO4+Ph9m39l1lAgA6t1vzZYh0xvwIFoHhxF27e1HD8uEJycgp5GVLGEXME/VxK1apVo3v37ly/fp2oqChatWrFxIkTGTFiBE5Opv1/3FHzsxRr5mfqPjWKYvpAx9euXaNMmTLs2bOHoKAgw/KRI0eyY8eOTEfAi4iI4OzZs9SsWZOYmBimT5/Ozp07OXHiBGXLls30OBMmTGDixIkZlq9du1ZuPRRPpAsXLrB582a2b9/OgwcPjN7z8fGhffv2tG/fnsKFC9spQnVKTYXz5304d84XNzcdzz2nglHPHpo4sRGHD+tbqZcu3Urx4vftHJH13Lt3jxkzZhgV8HXq1OG9996T/lNWlpiYSM+ePYmJicn2Z231YuRxWq2WKlWq0KNHDyZPnpzpOpm1jPj7+3Pr1i2L/uJotVrCw8Np06aNKgd5Unt+oO4cHzx4wDfffMP06dM5ffp0hvcbNGjAm2++yWuvvZavO4Wr+RyC4+Y3bpwT06bpZ17+9tsUunbN/QQ8jppjeqmpqUyZMoUpU6aQ9rVXtmxZVq9eTZMmTbLdNj/klxfWzC82NpZixYrlWIyYdZmmWLFiODs7Ex0dbbQ8OjqaUqVKmbQPV1dXAgMDiYyMzHIdd3d33N3dM93WGr8I1tqvo1B7fqCuHI8fP87KlSsJCwvjzp07Ru8VLFiQXr16MWjQINVNXKemc5gZR8uvQYNHz48edaFbt7zv09FyTM/V1ZXJkyfTvHlzevXqxc2bN7ly5QqtW7dm0qRJjBo1KsfLNo6cnyVYIz9T92dWB1Y3Nzfq1q3Ltm3bDMt0Oh3btm0zainJTmpqKseOHcMvbXIEIQR3795lwYIF1K9fnxo1ajBjxgyjQqRatWrMnz+fq1evsnjxYtUVIsL26tV79PzAAfvFYWtt2rTh8OHDhmHkU1NTGTNmDO3bt8/wH21hO2b3DQ8JCaF3797Uq1ePBg0aMGvWLBISEujbty8AwcHBlClThtDQUAAmTZpEo0aNqFChAvfu3ePzzz/n4sWLDBgwwLKZCJHPaLVatm7dypdffslPP/1kdGkS9MX/q6++So0aNQgJCcHNzc1OkT6Z4uJg/nz9l3bduqC27jj+/lCsmH4SwIMH9bMTPynzIpYpU4Zt27YxefJkJk2ahKIohIeHU6tWLcLCwmjfvr29Q3zimH1r7+uvv8706dMZP348tWvX5siRI2zevNlwu+6lS5e4fv26Yf27d+8ycOBAqlSpQseOHYmNjWXPnj1UrVrVclkIkU8oisKhQ4d4//33KVu2LC+88ALffvutUSFSp04d5s6dy7Vr1wgLC6NKlSoye64dHDyoH4ujTRsYM8be0VieRqMvskBfkKhhRmJzODs7M2HCBH7//XdDN4Po6Gg6dOjAsGHDuH9fvR16HVGu7pofOnQoQ4cOzfS97du3G72eOXMmM2fOzM1hhFCNyMhIvv76a7766qtMO6MWL16cXr160bdvX2rWrGlYrtZbCfOD9Jcu0l/SUJO6dWHLFv3zgwehXDn7xmMPrVq14ujRo/Tp08cwavGcOXPYtm0bX331FbVq1bJzhE8GlQzhI4TjuXz5Mt9//z3r1q3LdMoEd3d3XnzxRYKDg2nXrp2qO8blR09CMdK6NURHP7oU9aQqUaIEv/32G/Pnz2fEiBE8ePCAEydO0KBBAyZPnsy7775r7xBVT4oRISzo6tWrfP/993z77bfs2bMn03VatGhBz5496datm4wL4sDSihEPD1DrVeXnntM/hH6k8KFDh/Lcc8/Rq1cvjh49SnJyMh9++CE//PADb7zxhr1DVDUpRoTIo7QC5LvvvmP37t2ZrlOrVi169epF9+7d8ff3t3GEwlx378K5c/rngYHqGQZe5KxatWr8/fffjB07li+++AJFUdi3bx+HDx8mLi6OESNG4CK/EBYnP1EhcuG///5jw4YN/Pjjj1kWIDVq1KBbt25069aNypUr2zhCkRcHDz56Xr++/eIQ9uHu7s7nn3/OSy+9RN++fTlz5gxarZaPPvqIDRs2sGLFimynNBHmk2JECBMoisKRI0fYsGED69ev59ixY5muV61aNbp168Zrr71GlSpVbBylsJQnob9IenfvwqFDULYsVKpk72gcR+PGjTly5Ahjxoxh1qxZhlaSOnXqMH78eEaOHCl9vSxEihEhspCSksJff/3Fhg0b2LBhA5cuXcp0vbQCpFu3bnLLukrs3//oudqLkU2boGNH/fOPPoIpU+wbj6MpUKAA06ZNo2TJkqxYsYIzZ86QnJzM2LFj+eGHH1i5cqXccWMBUowIkc7NmzfZsmULGzduZPPmzdy9ezfT9Ro2bMhLL71E165dqST/lVSdtJaRQoXg2WftG4u1pb/a8CSNxGquypUrc+DAAaZMmcLnn3+OTqfj8OHD1KtXj1GjRvHRRx9RoEABe4eZb0kxIp5oOp2OAwcOsGnTJjZt2sS+ffvIbO5IFxcXWrVqRdeuXenSpQulS5e2Q7TCFrRa/S2vBw5A8eLg7GzviKzrSR6J1VweHh5MnTqVV155hb59+3LixAlSUlL45JNP+Prrr1m4cCFt2rSxd5j5khQj4olz69YtwsPD2bRpE5s3b+bmzZuZrufl5UWHDh3o2rUrHTp0wNfX17aBCrtwdYXly/XPdTr7xmILaSOxbtkCt2/DpUsQEGDvqBxb/fr1OXjwIFOmTGHq1KlotVrOnTtH27Zt6dGjBzNmzDB58lihJ8WIUD2tVktERARbtmxh69atHDx4MNPWD9DfAdOhQwc6duxI48aNpXPaEy6HSVxVo14945FYpRjJmbu7O5MmTaJ79+4MGjSIv/76C4Cvv/6aTZs2MXXqVAYOHJjjTMBCT4oRoTqKohAZGcnWrVvZunUrf/zxB/Hx8ZmuW6hQIVq3bk3Hjh1p3769jAEinkjpR189cABeftl+seQ3VatWZfv27YSFhTFixAju3LnDvXv3GDRoEF9++SWLFy+mevXq9g7T4UkxIlTh2rVr7Nq1iz///JNt27Zx4cKFLNetVasWbdq0oUOHDjRt2lRmwxVGYmLAx8feUdhW+mIk/RgrwjROTk7069ePzp07M2LECFatWgXAnj17CAwMZMSIEYwdOxZPT087R+q4pBgR+VJ0dDQ7duxg27Zt/Pbbb1y9ejXLdUuUKEGbNm1o164drVu3xs/Pz4aRivzk2jUoUwYqVoS+ffWz9j4JpBOrZRQvXpywsDB69+7NoEGDOHPmDCkpKYSGhrJ27VqmT5/OK6+8IrNwZ0KKEZEvXL9+nb/++osdO3bw559/curUqSzXdXNzo1mzZrRr1462bdtSo0YNuW4rTJJ2a+vZs/oWkieFdGK1rOeee46jR48ydepUQkNDSU5O5uLFi3Tr1o2WLVsye/Zso9m5hRQjwgEpisLFixf566+/2LlzJzt27ODs2bNZru/s7EyDBg1o1aoVrVq1IigoSO73F7ny99+Pnqt9sLPHpXVibd0aUlLsHU3+5+HhwYQJE+jevTvvvvsu4eHhAGzfvp3AwEDefPNNPvnkE4oWLWrnSB2DFCPC7nQ6HadOnWLXrl3s3LmTv/76i8uXL2e5vrOzM/Xq1aNly5Y0bdqU+Ph4XnnlFbnzReSJTgdff61/rtFAUJB947G1gQNhwAB46il7R6IulStXZsuWLfzyyy+EhIRw7tw5dDodixYt4ptvvuGTTz7hzTfffOIn33uysxd2kZiYyIEDB9i9eze7du1iz5493Lt3L8v1XV1dadCgAc2bN6dZs2Y0adIEb29vQH/b7saNG20UuVCzP/+E8+f1z9u00fcdeZLIZRnr0Wg0vPjii7Rr145Zs2bxySefEB8fz927dxkyZAhLlixhzpw5NG/e3N6h2o0UI8KqFEXhwoULREREGB5Hjx4lJZt2YE9PT4KCggzFR6NGjeSyi7C6ZcsePR8wwH5xCPVyd3fnww8/5I033mDUqFGsXr0agKNHj9KiRQtee+01Pv30U5555hk7R2p7UowIi4qPj+fgwYPs3buXiIgI9u7dS3R0dLbbFC9enCZNmtCkSROaNWtGnTp15JKLsKlbt+DHH/XPixWDLl3sG4+9pabCtm3QqhU84VcPrKJ06dJ8+eWXvPXWW7z77rscOnQIgG+//Zb169czePBgxo4dS/Hixe0cqe3Ir5nIteTkZI4dO8a+ffvYv38/+/bt49SpU+iyGUNbo9FQtWpVgoKCaNy4MU2aNKFixYpyq5uwqzVrIDlZ/7x3b3iSh5755hsYMQIuX4ZffoEXXrB3ROrVpEkT9u3bx/Llyxk3bhw3btxAq9UyZ84cVq5cyahRo3jvvfeeiPFJpBgRJklNTeXMmTMcPHjQUHgcPnyYpKSkbLfz8fGhUaNGNGrUiKCgIBo2bChzvAiHoiiwdOmj1/372y8WR1CwoL4QAf0cPVKMWJezszNvvvkmPXr04PPPP+eLL74gMTGRuLg4xowZw/z585k8eTK9e/fGWcWzNkoxIjLQ6XScOXOGAwcOcPDgQQ4ePMjhw4ezHFI9jYuLCzVr1qR+/fo0aNCAoKAgKlWqJGN8CId2/z40a6b/Aq5VC6pUsXdE9tW+Pfj5wfXr8OuvEB0NJUvaOyr18/LyYtKkSQwePJgJEyawbNkydDod165do3///sycOZNp06bRoUMHVbYkSzHyhNNqtZw6dYojR45w5MgRDh48yKFDh3IsPAAqVqxIgwYNqF+/Pg0bNqRWrVrS0VTkO56esGgRfPGF/gv4Sefior9UNXWqfryRNWtg+HB7R/Xk8PPzY/Hixbz33nuMGjWKn3/+GYDjx4/TqVMnWrZsSWhoKI0aNbJzpJYlxcgTJDY2lqNHjxoKjyNHjnD8+HGS0y6WZ6NcuXLUq1ePunXrUq9ePerVq0eRIkVsELUQtlGwIFSoYO8oHEO/fvpiBPSXakJCZHh4W6tSpQo//fQTO3fuZMSIEezbtw/QD5oWFBREly5d+OSTT1QzCZ8UIyqk0+m4cOECBw8e5Mcff2TlypUcO3aM//77z6Tty5UrR926dY0eT1KvbiGedBUr6i9d/fUXnDqlH5lWZf8RzzeaN2/O3r17+e677xgzZgyRkZEA/PTTT/z888/873//Y8KECTz99NN2jjRvpBjJ52JiYjh+/Dj//POP4XHs2DHi4uJy3Faj0VCpUiUCAwOpXbu24VGiRAkbRC6EfSmK/hJE584gfaoz6tdPX4wAzJsnxYg9aTQaXnvtNV566SVWrFjBpEmTuHbtGoqisHr1ar7++mv69evH2LFj8ff3t3e4uSLFSD5x//59/v33X44fP254HDt2LNth09Pz9PSkRo0aRoVH9erVKViwoJUjF8IxTZgAkyaBhwd89BGMG2fviBxLt27w7rsQFwdffQXNm8Obb9o7qiebq6srb731FsHBwcybN4/Q0FDu3r1LSkoKS5YsISwsjEGDBjF69GhKlSpl73DNIsWIg3nw4AGnT5/mxIkTnDx5khMnTnDixAnDfAameOqpp6hZsybVq1cnJSWF3r17U7lyZbmrRYiHVq3SFyIASUn6GWuFsYIFYcYM/Zw1ACtW6G97VvHdpflGgQIFGDFiBG+++SYzZ85kxowZxMXFkZyczJw5c1i6dClDhgxh5MiR+eYSuxQjdpKQkMC///7LyZMnOXnyJKdOneLkyZNmFR3e3t7UqFGD6tWrU6NGDWrVqkWNGjXw8fEBHs3bUrFiRSlEhHjozz8ffcGC/i6ajh3tF48jGzBA32fk/HlYvVoKEUfj4+PDhAkTeOedd/j888+ZM2cO9+/f5/79+0yfPp1Fixbx7rvv8sEHH1C4cGF7h5stKUasSFEUbt68yb///mt4pBUdly5dMnk/BQoUoGrVqlSrVo3q1atTrVo1atSoQdmyZVV5v7kQ1nLqFLz0Emi1+tdDhsB779k1JIf32Wf6O2nk/zOOq2jRokydOpX333+fqVOnsnDhQpKSkoiPj+fTTz9l3rx5DBs2jPfee89h74KUYsQCkpOTOXfuHKdPn+bff/81+vfu3bsm76dAgQJUrlyZatWqGR5Vq1blqaeeUvXIe0LYwvHj+s6qMTH61506waxZcstqTjL703Pzpr4vST6/gUN1SpYsycyZM/nggw/49NNPWbp0KVqtltjYWCZPnsysWbN45513CAkJoWjRovYO10iuipH58+fz+eefExUVRa1atZg7dy4NGjTIcv3vvvuOcePGceHCBSpWrMi0adPomM/aRVNTU7l8+TJnz57lzJkzhsfZs2c5f/68yZdWQN+0VqVKFcOjWrVqVKlShYCAALmcIoSFXb+u73j566+PlgUGwrp1Mglcbly4oC/qTp2CHj1g1Ch49ll7RyXSK1OmDPPnz2fkyJFMmTKFlStXkpKSQlxcHJ9++ilz5sxh6NChhISEOEyfErM/it988w0hISEsWrSIhg0bMmvWLNq1a8fp06czvSV0z5499OjRg9DQUF544QXWrl1L165dOXTokMMN1qLT6bh69SqRkZGcPXuWs2fPGp5HRkbmOA/L4/z9/alSpQqVKlWicuXKVKpUiapVq1KqVCm5vCKEFSiKvuUj/a26hQvDgQOPXpcvry9MChWyeXj5XlISdOgA//6rf71mTdrt0c6ULRuAr6+GWrXkVmlHERAQwJIlSxgzZgyhoaGsWLECrVZLfHw8U6dOZe7cubz99tsMGzbM3qGiURRFMWeDhg0bUr9+febNmwfov8D9/f155513GDVqVIb1X3/9dRISEvg13X9LGjVqRO3atVm0aJFJx4yNjcXHx4eYmBi8vb3NCTdL0dHRrF+/nuXL40lIKM7t27e5c+cOKSnaLLbYD/yZ7rUGGIGbmzvFihWjePFiFC9ewvBvsWLFcHs49edrr+n/AKY5dw6+/z7nGDUaGDnSeNnvv8PBgzlv+/TT0LWrvgNrx44dcXV1ZelSuHMn522ffx7q1Xv0OiZGP1y2KQYMgPStf0eOwJYtOW/n7Q2DBxsvW78ezpzJfrvU1FSSkvYxdmw9XF1dDctnzHjULyA7XbtCpUqPXl+5or+N0RTvv288u+tff8GePTlvV6YM/O9/xsu+/DLzochTU1M5ffpfKlWqjLOzM02bQpMmj95/8ABmzzYt3jfegNKlH70+dQoejjSdLXf3jP0qNm6EY8dy3rZyZejSxXjZvHmQfraBlJRUTp8+TaVKlQyXIxVFP0FbzZqP1rt+HRYs0A9RnpKi/2KMi3v0iImByEj9eb93z/jyy9Sp+m1DQvS/o7YsRNI6kqd9DvO7u3dh/nz9Ja7btzNfp0wZ/TglaXcspRk3Tn/uXFz0j7SGYI3G+PHyy8afy8uXc/5cpp3vkBBI/2PeuRMiInLOq2xZ6NXLeNmqVRAVlfFz+LhmzaBx40ev79+HOXNyPiZAcLB+LqA0J0/qZ0vOiYcHPF5D/Pab/nJkVu7du8v27Ts4cGA1KSk/GpZ7enry3HOdWLRoBmXLljUtcBOZ/P2tmCEpKUlxdnZW1q9fb7Q8ODhYefHFFzPdxt/fX5k5c6bRsvHjxys1a9bM8jgPHjxQYmJiDI/Lly8rgHLr1i0lOTnZIo+///5bARTYquj/9GX/KFJkhfLiiy8qw4cPVxYtWqRs3brNpO1AUTZu1Bod+6eftCZt5+SkyxD30KEpJm3bvn2qkpCQoGzYsEFJSEhQkpOTlWef1Zm07RdfpBgd8+zZZJNz/ecf43gXLDAt13LlMubatWuqSdu2aXPekGPao2BB03Jdt8743Gzfblq8oCh37hjHO3q0aecmKCg1Q6716pmW68cfG5+b6GjTz83u3ca5rl5tWq6+vhnPTXCwafF265YxVz8/087N0qXG8e7fb3qu588bHzMmJllJSLDM3w5zH49/DtXyuHMnWfn88xSldOnMz2e/fhnPvaenaef+22+Nz/0ff5j+ubx3z/iYo0aZ9rls0iRjvHXrmvZ7PmGC8ecyKsr039U9e4xz/fJL03ItXDjj5/KNN0yL94UX4pS3335bcXd3V/TfgyhQVKlUqZKSlJRk0d+TW7duKYASExOTbX1h1mWaW7dukZqaSsnHpnAsWbIk/6a12z0mKioq0/WjoqKyPE5oaCgTJ07MsHzr1q14enqaE3KWEhMTzVq/WbNm9O37qBeyKSOcptm3bx/JyTcNrw8cKAEEmbTtxo0bjV5fuFAdeCbH7W7cuEl4+F4AwsPDAUhIaAV45bjtyZMn2bjx0dDxN24UANqaFO/OnTv5779H/+09diwAqJ3jdvfv32fjxnCjZVFR9YHSmW/wmLQc06SmdsKUq5CHDh3Cw+NRk8SpU0WAZiYdc8uWLRQokGp4HRlZGaiU9QYP3b17l40bdxkti4lpDuR8692ZM2fYuPFRc1F8vCtgWv+r3bt3c/PmPcPrI0fKAPWyXD+N/n/2m4yWXbkSCJTLcdvr16+xcaNxU96DB+0Ajxy3/eefY2zc+Oius0uXvIBW2W5TpMh9/P3j2LjxGGXK5DzZoy09/juqBhUrwuzZTpw+XZhLl7y5dMmLS5e8uHzZm1u3LrFx4wmj9ZOTXwBy7ox/6NAh3NwefS5PnDDvc+nu/uhzee6caZ/LO3cs97mMizPvc3njxj3Da1t8LmNi7tG2bVvq1avHjz/+yNatW0lOhlatWrFp06YctzeHqd+1Zl2muXbtGmXKlGHPnj0EBT36Mh05ciQ7duzg77//zrCNm5sbq1atokePHoZlCxYsYOLEiURHR2d6nKSkJKP+GbGxsfj7+3Pr1i2LXaYBWLZsGQcOuFO7dhtKliyZbefRZ55RSN/FRVHg559N6/cRFKSQvjtNVBT8/bdp23bpYnx6/vkHzp/PeduSJaFu3WTCw8Np06YNrq6ubNumwYTJeKleXeGZdPVOYiKEh5sWb6tWCl7p6p0LF+Do0Zy3LVAA2rY1znXfPk2Os6impqZy7VoEb73VwKgJ/NdfNaSmZrPhQ/XrK0aXLm7fhl27TMu1UyfFqAPkqVNw5kzO2xYpAs2aGee6Y4eGe/cyrpuamsqRI0eoXbs2zs7OVKqkULnyo/eTk2HTJtPibd5cIf1QA1euwMGDOW/r6godOxrHe/gwXLqU87ZlykC9esbbbtmiIX33q8dzTFO7tkK5dH9X4+P1nxsXF/0dHm5uUKiQ/vfNy0t/6cURO6RqtVqjz6EaZZajomS8U2n3bo3hMptWCzpdxv+3g/5zmf7Sxa1b+m1N8fjn8uRJOHs2522LFoWmTTP/XGb1O5qmcmXF6LJSfvtcXrp0iQ8/HM/SpfMoZOFrmLGxsRQrVix/XqZ5XExMjEnNPOZKTk5WNmzYoCQnJ1t0v45C7fkpivpzVHt+iqL+HNWen6KoP0fJL/dM/f426z5SNzc36taty7Zt2wzLdDod27ZtM2opSS8oKMhofdA3V2a1vhBCCCGeLGY3aoaEhNC7d2/q1atHgwYNmDVrFgkJCfTt2xeA4OBgypQpQ2hoKADDhg2jRYsWfPHFF3Tq1Il169Zx4MABlixZYtlMhBBCCJEvmV2MvP7669y8eZPx48cTFRVF7dq12bx5s6GT6qVLl4z6XjRu3Ji1a9cyduxYPvroIypWrMiGDRscbowRIYQQQthHrrp7DR06lKFDh2b63vbt2zMs69atG926dcvNoYQQQgihcjL2uBBCCCHsSooRIYQQQtiVFCNCCCGEsCspRoQQQghhV1KMCCGEEMKupBgRQgghhF1JMSKEEEIIu5JiRAghhBB2JcWIEEIIIezKASfczkh5OK90bGysRfer1WpJTEwkNjZWlVN7qz0/UH+Oas8P1J+j2vMD9eco+eVe2vd22vd4VvJFMRIXFweAv7+/nSMRQgghhLni4uLw8fHJ8n2NklO54gB0Oh3Xrl3Dy8sLjUZjsf3Gxsbi7+/P5cuX8fb2tth+HYXa8wP156j2/ED9Oao9P1B/jpJf7imKQlxcHKVLlzaaRPdx+aJlxMnJibJly1pt/97e3qr8BUuj9vxA/TmqPT9Qf45qzw/Un6PklzvZtYikkQ6sQgghhLArKUaEEEIIYVdPdDHi7u7Oxx9/jLu7u71DsQq15wfqz1Ht+YH6c1R7fqD+HCU/68sXHViFEEIIoV5PdMuIEEIIIexPihEhhBBC2JUUI0IIIYSwKylGhBBCCGFXqi5GpkyZQuPGjfH09MTX19ekbRRFYfz48fj5+VGgQAFat27N2bNnjda5c+cOvXr1wtvbG19fX/r37098fLwVMsiZubFcuHABjUaT6eO7774zrJfZ++vWrbNFSkZy87Nu2bJlhtgHDRpktM6lS5fo1KkTnp6elChRghEjRpCSkmLNVLJkbo537tzhnXfeoVKlShQoUIBy5crx7rvvEhMTY7Sevc7h/Pnzeeqpp/Dw8KBhw4bs27cv2/W/++47KleujIeHBzVq1GDjxo1G75vymbQ1c3JcunQpzZo1o3DhwhQuXJjWrVtnWL9Pnz4ZzlX79u2tnUaWzMkvLCwsQ+weHh5G6+T3c5jZ3xSNRkOnTp0M6zjSOdy5cyedO3emdOnSaDQaNmzYkOM227dvp06dOri7u1OhQgXCwsIyrGPuZ9ssioqNHz9emTFjhhISEqL4+PiYtM3UqVMVHx8fZcOGDcrRo0eVF198USlfvrxy//59wzrt27dXatWqpezdu1f566+/lAoVKig9evSwUhbZMzeWlJQU5fr160aPiRMnKoUKFVLi4uIM6wHKypUrjdZL/zOwldz8rFu0aKEMHDjQKPaYmBjD+ykpKUr16tWV1q1bK4cPH1Y2btyoFCtWTBk9erS108mUuTkeO3ZMefnll5Wff/5ZiYyMVLZt26ZUrFhReeWVV4zWs8c5XLduneLm5qasWLFCOXHihDJw4EDF19dXiY6OznT93bt3K87Ozspnn32mnDx5Uhk7dqzi6uqqHDt2zLCOKZ9JWzI3x549eyrz589XDh8+rJw6dUrp06eP4uPjo1y5csWwTu/evZX27dsbnas7d+7YKiUj5ua3cuVKxdvb2yj2qKgoo3Xy+zm8ffu2UX7Hjx9XnJ2dlZUrVxrWcaRzuHHjRmXMmDHKjz/+qADK+vXrs13/v//+Uzw9PZWQkBDl5MmTyty5cxVnZ2dl8+bNhnXM/ZmZS9XFSJqVK1eaVIzodDqlVKlSyueff25Ydu/ePcXd3V35+uuvFUVRlJMnTyqAsn//fsM6mzZtUjQajXL16lWLx54dS8VSu3ZtpV+/fkbLTPkFtrbc5teiRQtl2LBhWb6/ceNGxcnJyegP5sKFCxVvb28lKSnJIrGbylLn8Ntvv1Xc3NwUrVZrWGaPc9igQQNlyJAhhtepqalK6dKlldDQ0EzXf+2115ROnToZLWvYsKHy1ltvKYpi2mfS1szN8XEpKSmKl5eXsmrVKsOy3r17K126dLF0qLlibn45/X1V4zmcOXOm4uXlpcTHxxuWOdI5TM+UvwMjR45UqlWrZrTs9ddfV9q1a2d4ndefWU5UfZnGXOfPnycqKorWrVsblvn4+NCwYUMiIiIAiIiIwNfXl3r16hnWad26NU5OTvz99982jdcSsRw8eJAjR47Qv3//DO8NGTKEYsWK0aBBA1asWJHjFNCWlpf8vvrqK4oVK0b16tUZPXo0iYmJRvutUaMGJUuWNCxr164dsbGxnDhxwvKJZMNSv08xMTF4e3vj4mI83ZQtz2FycjIHDx40+vw4OTnRunVrw+fncREREUbrg/5cpK1vymfSlnKT4+MSExPRarUUKVLEaPn27dspUaIElSpVYvDgwdy+fduisZsit/nFx8cTEBCAv78/Xbp0MfocqfEcLl++nO7du1OwYEGj5Y5wDnMjp8+hJX5mOckXE+XZSlRUFIDRl1Ta67T3oqKiKFGihNH7Li4uFClSxLCOrVgiluXLl1OlShUaN25stHzSpEm0atUKT09Ptm7dyttvv018fDzvvvuuxeLPSW7z69mzJwEBAZQuXZp//vmHDz/8kNOnT/Pjjz8a9pvZOU57z5YscQ5v3brF5MmTefPNN42W2/oc3rp1i9TU1Ex/tv/++2+m22R1LtJ/3tKWZbWOLeUmx8d9+OGHlC5d2ugPe/v27Xn55ZcpX748586d46OPPqJDhw5ERETg7Oxs0Ryyk5v8KlWqxIoVK6hZsyYxMTFMnz6dxo0bc+LECcqWLau6c7hv3z6OHz/O8uXLjZY7yjnMjaw+h7Gxsdy/f5+7d+/m+fc+J/muGBk1ahTTpk3Ldp1Tp05RuXJlG0VkeabmmFf3799n7dq1jBs3LsN76ZcFBgaSkJDA559/bpEvMmvnl/5LuUaNGvj5+fH8889z7tw5nnnmmVzv1xy2OoexsbF06tSJqlWrMmHCBKP3rHkORe5MnTqVdevWsX37dqNOnt27dzc8r1GjBjVr1uSZZ55h+/btPP/88/YI1WRBQUEEBQUZXjdu3JgqVaqwePFiJk+ebMfIrGP58uXUqFGDBg0aGC3Pz+fQEeS7YmT48OH06dMn23WefvrpXO27VKlSAERHR+Pn52dYHh0dTe3atQ3r3Lhxw2i7lJQU7ty5Y9g+r0zNMa+xfP/99yQmJhIcHJzjug0bNmTy5MkkJSXlef4CW+WXpmHDhgBERkbyzDPPUKpUqQy9wKOjowHy1TmMi4ujffv2eHl5sX79elxdXbNd35LnMDPFihXD2dnZ8LNMEx0dnWUupUqVynZ9Uz6TtpSbHNNMnz6dqVOn8vvvv1OzZs1s13366acpVqwYkZGRNv0iy0t+aVxdXQkMDCQyMhJQ1zlMSEhg3bp1TJo0Kcfj2Osc5kZWn0Nvb28KFCiAs7Nznn8vcmSRnicOztwOrNOnTzcsi4mJybQD64EDBwzrbNmyxa4dWHMbS4sWLTLcgZGVTz75RClcuHCuY80NS/2sd+3apQDK0aNHFUV51IE1fS/wxYsXK97e3sqDBw8sl4AJcptjTEyM0qhRI6VFixZKQkKCSceyxTls0KCBMnToUMPr1NRUpUyZMtl2YH3hhReMlgUFBWXowJrdZ9LWzM1RURRl2rRpire3txIREWHSMS5fvqxoNBrlp59+ynO85spNfumlpKQolSpVUt5//31FUdRzDhVF/13i7u6u3Lp1K8dj2PMcpoeJHVirV69utKxHjx4ZOrDm5fcixzgtshcHdfHiReXw4cOGW1cPHz6sHD582OgW1kqVKik//vij4fXUqVMVX19f5aefflL++ecfpUuXLpne2hsYGKj8/fffyq5du5SKFSva9dbe7GK5cuWKUqlSJeXvv/822u7s2bOKRqNRNm3alGGfP//8s7J06VLl2LFjytmzZ5UFCxYonp6eyvjx462ez+PMzS8yMlKZNGmScuDAAeX8+fPKTz/9pDz99NNK8+bNDduk3drbtm1b5ciRI8rmzZuV4sWL2/XWXnNyjImJURo2bKjUqFFDiYyMNLqVMCUlRVEU+53DdevWKe7u7kpYWJhy8uRJ5c0331R8fX0Ndy698cYbyqhRowzr7969W3FxcVGmT5+unDp1Svn4448zvbU3p8+kLZmb49SpUxU3Nzfl+++/NzpXaX+H4uLilA8++ECJiIhQzp8/r/z+++9KnTp1lIoVK9q8OM5NfhMnTlS2bNminDt3Tjl48KDSvXt3xcPDQzlx4oRhnfx+DtM0bdpUef311zMsd7RzGBcXZ/i+A5QZM2Yohw8fVi5evKgoiqKMGjVKeeONNwzrp93aO2LECOXUqVPK/PnzM721N7ufWV6puhjp3bu3AmR4/Pnnn4Z1eDgWQxqdTqeMGzdOKVmypOLu7q48//zzyunTp432e/v2baVHjx5KoUKFFG9vb6Vv375GBY4t5RTL+fPnM+SsKIoyevRoxd/fX0lNTc2wz02bNim1a9dWChUqpBQsWFCpVauWsmjRokzXtTZz87t06ZLSvHlzpUiRIoq7u7tSoUIFZcSIEUbjjCiKoly4cEHp0KGDUqBAAaVYsWLK8OHDjW6LtSVzc/zzzz8z/b0GlPPnzyuKYt9zOHfuXKVcuXKKm5ub0qBBA2Xv3r2G91q0aKH07t3baP1vv/1WefbZZxU3NzelWrVqym+//Wb0vimfSVszJ8eAgIBMz9XHH3+sKIqiJCYmKm3btlWKFy+uuLq6KgEBAcrAgQMt9kc+N8zJ77333jOsW7JkSaVjx47KoUOHjPaX38+hoijKv//+qwDK1q1bM+zL0c5hVn8j0nLq3bu30qJFiwzb1K5dW3Fzc1Oefvppo+/FNNn9zPJKoyg2vl9TCCGEECIdGWdECCGEEHYlxYgQQggh7EqKESGEEELYlRQjQgghhLArKUaEEEIIYVdSjAghhBDCrqQYEUIIIYRdSTEihBBCCLuSYkQIIYQQdiXFiBBCCCHsSooRIYQQQtiVFCNCCJv78ssvKVq0KElJSUbLu3btyhtvvGGnqIQQ9iLFiBDC5rp160Zqaio///yzYdmNGzf47bff6Nevnx0jE0LYgxQjQgibK1CgAD179mTlypWGZWvWrKFcuXK0bNnSfoEJIexCihEhhF0MHDiQrVu3cvXqVQDCwsLo06cPGo3GzpEJIWxNoyiKYu8ghBBPprp16/Lqq6/Stm1bGjRowIULF/D397d3WEIIG3OxdwBCiCfXgAEDmDVrFlevXqV169ZSiAjxhJKWESGE3cTExFC6dGlSUlL48ssvef311+0dkhDCDqTPiBDCbnx8fHjllVcoVKgQXbt2tXc4Qgg7kWJECGFXV69epVevXri7u9s7FCGEnchlGiGEXdy9e5ft27fz6quvcvLkSSpVqmTvkIQQdiIdWIUQdhEYGMjdu3eZNm2aFCJCPOGkZUQIIYQQdiV9RoQQQghhV1KMCCGEEMKupBgRQgghhF1JMSKEEEIIu5JiRAghhBB2JcWIEEIIIexKihEhhBBC2JUUI0IIIYSwq/8D2FIJS+BTgiYAAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "fig = plt.figure()\n", @@ -1289,7 +434,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1332,7 +477,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/chapter1/membrane_code.py b/chapter1/membrane_code.py index 7a217a82..770bf66a 100644 --- a/chapter1/membrane_code.py +++ b/chapter1/membrane_code.py @@ -27,7 +27,7 @@ # # ## Creating the mesh # -# To create the computational geometry, we use the python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it. +# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it. import gmsh gmsh.initialize() @@ -171,8 +171,8 @@ def on_boundary(x): bb_tree = geometry.bb_tree(domain, domain.topology.dim) # Now we can compute which cells the bounding box tree collides with using `dolfinx.geometry.compute_collisions_points`. This function returns a list of cells whose bounding box collide for each input point. As different points might have different number of cells, the data is stored in `dolfinx.cpp.graph.AdjacencyList_int32`, where one can access the cells for the `i`th point by calling `links(i)`. -# However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell, we check that the actual cell collides with cell -# using `dolfinx.geometry.select_colliding_cells`, who measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries). +# However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell, we check that the actual cell collides with cell +# using `dolfinx.geometry.select_colliding_cells`, which measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries). # This function also returns an adjacency-list, as the point might align with a facet, edge or vertex that is shared between multiple cells in the mesh. # # Finally, we would like the code below to run in parallel, when the mesh is distributed over multiple processors. In that case, it is not guaranteed that every point in `points` is on each processor. Therefore we create a subset `points_on_proc` only containing the points found on the current processor. @@ -188,7 +188,8 @@ def on_boundary(x): points_on_proc.append(point) cells.append(colliding_cells.links(i)[0]) -# We now got a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points. +# We now have a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points. +# points_on_proc = np.array(points_on_proc, dtype=np.float64) u_values = uh.eval(points_on_proc, cells) diff --git a/chapter2/diffusion_code.ipynb b/chapter2/diffusion_code.ipynb index c468c004..a55acfa2 100644 --- a/chapter2/diffusion_code.ipynb +++ b/chapter2/diffusion_code.ipynb @@ -51,7 +51,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that we have used a much higher resolution that before to better resolve features of the solution.\n", + "Note that we have used a much higher resolution than before to better resolve features of the solution.\n", "We also easily update the intial and boundary conditions. Instead of using a class to define the initial condition, we simply use a function" ] }, @@ -82,7 +82,7 @@ "metadata": {}, "source": [ "## Time-dependent output\n", - "To visualize the solution in an external program such as Paraview, we create a an `XDMFFile` which we can store multiple solutions in. The main advantage with an XDMFFile, is that we only need to store the mesh once, and can append multiple solutions to the same grid, reducing the storage space.\n", + "To visualize the solution in an external program such as Paraview, we create a an `XDMFFile` which we can store multiple solutions in. The main advantage with an XDMFFile is that we only need to store the mesh once and that we can append multiple solutions to the same grid, reducing the storage space.\n", "The first argument to the XDMFFile is which communicator should be used to store the data. As we would like one output, independent of the number of processors, we use the `COMM_WORLD`. The second argument is the file name of the output file, while the third argument is the state of the file,\n", "this could be read (`\"r\"`), write (`\"w\"`) or append (`\"a\"`)." ] @@ -164,7 +164,7 @@ "metadata": {}, "source": [ "## Using petsc4py to create a linear solver\n", - "As we have already assembled `a` into the matrix `A`, we can no longer use the `dolfinx.fem.petsc.LinearProblem` class to solve the problem. Therefore, we create a linear algebra solver using PETSc, and assign the matrix `A` to the solver, and choose the solution strategy." + "As we have already assembled `a` into the matrix `A`, we can no longer use the `dolfinx.fem.petsc.LinearProblem` class to solve the problem. Therefore, we create a linear algebra solver using PETSc, assign the matrix `A` to the solver, and choose the solution strategy." ] }, { @@ -224,8 +224,8 @@ "## Updating the solution and right hand side per time step\n", "To be able to solve the variation problem at each time step, we have to assemble the right hand side and apply the boundary condition before calling\n", "`solver.solve(b, uh.vector)`. We start by resetting the values in `b` as we are reusing the vector at every time step.\n", - "The next step is to assemble the vector, calling `dolfinx.fem.petsc.assemble_vector(b, L)` which means that we are assemble the linear for `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side.\n", - "This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ if the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions.\n", + "The next step is to assemble the vector calling `dolfinx.fem.petsc.assemble_vector(b, L)`, which means that we are assembling the linear form `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side.\n", + "This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ in the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions.\n", "When we have applied the boundary condition, we can solve the linear system and update values that are potentially shared between processors.\n", "Finally, before moving to the next time step, we update the solution at the previous time step to the solution at this time step." ] @@ -283,7 +283,7 @@ "\n", "Then, we add a time-annotation to the figure, pressing: `Sources->Alphabetical->Annotate Time` and `Apply` in the properties panel. It Is also a good idea to select an output resolution, by pressing `View->Preview->1280 x 720 (HD)`.\n", "\n", - "Then finally, click `File->Save Animation`, and save the animation to the desired format, such as `avi`, `ogv` or a sequence of `png`s. Make sure to set the framerate to something, sensible, in the range of $5-10$ frames per second." + "Then finally, click `File->Save Animation`, and save the animation to the desired format, such as `avi`, `ogv` or a sequence of `png`s. Make sure to set the frame rate to something sensible, in the range of $5-10$ frames per second." ] } ], diff --git a/chapter2/diffusion_code.py b/chapter2/diffusion_code.py index 70ce040c..7315cf29 100644 --- a/chapter2/diffusion_code.py +++ b/chapter2/diffusion_code.py @@ -52,7 +52,7 @@ # - -# Note that we have used a much higher resolution that before to better resolve features of the solution. +# Note that we have used a much higher resolution than before to better resolve features of the solution. # We also easily update the intial and boundary conditions. Instead of using a class to define the initial condition, we simply use a function # + @@ -73,7 +73,7 @@ def initial_condition(x, a=5): # - # ## Time-dependent output -# To visualize the solution in an external program such as Paraview, we create a an `XDMFFile` which we can store multiple solutions in. The main advantage with an XDMFFile, is that we only need to store the mesh once, and can append multiple solutions to the same grid, reducing the storage space. +# To visualize the solution in an external program such as Paraview, we create a an `XDMFFile` which we can store multiple solutions in. The main advantage with an XDMFFile is that we only need to store the mesh once and that we can append multiple solutions to the same grid, reducing the storage space. # The first argument to the XDMFFile is which communicator should be used to store the data. As we would like one output, independent of the number of processors, we use the `COMM_WORLD`. The second argument is the file name of the output file, while the third argument is the state of the file, # this could be read (`"r"`), write (`"w"`) or append (`"a"`). @@ -109,7 +109,7 @@ def initial_condition(x, a=5): b = create_vector(linear_form) # ## Using petsc4py to create a linear solver -# As we have already assembled `a` into the matrix `A`, we can no longer use the `dolfinx.fem.petsc.LinearProblem` class to solve the problem. Therefore, we create a linear algebra solver using PETSc, and assign the matrix `A` to the solver, and choose the solution strategy. +# As we have already assembled `a` into the matrix `A`, we can no longer use the `dolfinx.fem.petsc.LinearProblem` class to solve the problem. Therefore, we create a linear algebra solver using PETSc, assign the matrix `A` to the solver, and choose the solution strategy. solver = PETSc.KSP().create(domain.comm) solver.setOperators(A) @@ -143,8 +143,8 @@ def initial_condition(x, a=5): # ## Updating the solution and right hand side per time step # To be able to solve the variation problem at each time step, we have to assemble the right hand side and apply the boundary condition before calling # `solver.solve(b, uh.vector)`. We start by resetting the values in `b` as we are reusing the vector at every time step. -# The next step is to assemble the vector, calling `dolfinx.fem.petsc.assemble_vector(b, L)` which means that we are assemble the linear for `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side. -# This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ if the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions. +# The next step is to assemble the vector calling `dolfinx.fem.petsc.assemble_vector(b, L)`, which means that we are assembling the linear form `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side. +# This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ in the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions. # When we have applied the boundary condition, we can solve the linear system and update values that are potentially shared between processors. # Finally, before moving to the next time step, we update the solution at the previous time step to the solution at this time step. @@ -185,4 +185,4 @@ def initial_condition(x, a=5): # # Then, we add a time-annotation to the figure, pressing: `Sources->Alphabetical->Annotate Time` and `Apply` in the properties panel. It Is also a good idea to select an output resolution, by pressing `View->Preview->1280 x 720 (HD)`. # -# Then finally, click `File->Save Animation`, and save the animation to the desired format, such as `avi`, `ogv` or a sequence of `png`s. Make sure to set the framerate to something, sensible, in the range of $5-10$ frames per second. +# Then finally, click `File->Save Animation`, and save the animation to the desired format, such as `avi`, `ogv` or a sequence of `png`s. Make sure to set the frame rate to something sensible, in the range of $5-10$ frames per second. diff --git a/chapter2/heat_code.ipynb b/chapter2/heat_code.ipynb index bdd129a6..f59c9dbe 100644 --- a/chapter2/heat_code.ipynb +++ b/chapter2/heat_code.ipynb @@ -212,7 +212,7 @@ "metadata": {}, "source": [ "## Solving the time-dependent problem\n", - "With these structures in place, we crete our time-stepping loop.\n", + "With these structures in place, we create our time-stepping loop.\n", "In this loop, we first update the Dirichlet boundary condition by interpolating the updated\n", "expression `u_exact` into `V`. The next step is to re-assemble the vector `b`, with the update `u_n`.\n", "Then, we need to apply the boundary condition to this vector. We do this by using the lifting operation,\n", diff --git a/chapter2/heat_code.py b/chapter2/heat_code.py index 2586ae09..9603da9f 100644 --- a/chapter2/heat_code.py +++ b/chapter2/heat_code.py @@ -111,7 +111,7 @@ def __call__(self, x): solver.getPC().setType(PETSc.PC.Type.LU) # ## Solving the time-dependent problem -# With these structures in place, we crete our time-stepping loop. +# With these structures in place, we create our time-stepping loop. # In this loop, we first update the Dirichlet boundary condition by interpolating the updated # expression `u_exact` into `V`. The next step is to re-assemble the vector `b`, with the update `u_n`. # Then, we need to apply the boundary condition to this vector. We do this by using the lifting operation, diff --git a/chapter2/heat_equation.md b/chapter2/heat_equation.md index f66e03cc..d36aac28 100644 --- a/chapter2/heat_equation.md +++ b/chapter2/heat_equation.md @@ -9,7 +9,7 @@ As a first extension of the Poisson problem from the previous chapter, we consid The model problem for the time-dependent PDE reads \begin{align} \frac{\partial u}{\partial t}&=\nabla^2 u + f && \text{in } \Omega \times (0, T],\\ - u &= u_D && \text{n } \partial\Omega \times (0,T],\\ + u &= u_D && \text{on } \partial\Omega \times (0,T],\\ u &= u_0 && \text{at } t=0. \end{align} @@ -38,7 +38,7 @@ We reorder the equation such that the left-hand side contains the terms with onl \end{align} Given $u_0$, we can solve for $u^0, u^1, u^2$ and so on. -We then in turn use the finite element method. This means that we have to turn the equation into its weak formulation. We multiply by the test-function of $v\in \hat{V}$ and integrate second-order derivatives by parts. we now introduce the symbol $u$ for $u^{n+1}$ and we write the resulting weak formulation as +We then in turn use the finite element method. This means that we have to turn the equation into its weak formulation. We multiply by the test-function of $v\in \hat{V}$ and integrate second-order derivatives by parts. We now introduce the symbol $u$ for $u^{n+1}$ and we write the resulting weak formulation as \begin{align} a(u,v)&=L_{n+1}(v), @@ -64,4 +64,4 @@ When solving this variational problem $u^0$ becomes the $L^2$-projection of the The alternative is to construct $u^0$ by just interpolating the initial value $u_0$. We covered how to use interpolation in DOLFINx in the {doc}`membrane chapter <../chapter1/membrane_code>`. -We can use DOLFINx to either project or interpolate the initial condition. The most common choice is to use an projection, which computes an approximation to $u_0$. However, in some applications where we want to verify the code by reproducing exact solutions, one must use interpolate. In this chapter, we will use such a problem. +We can use DOLFINx to either project or interpolate the initial condition. The most common choice is to use a projection, which computes an approximation to $u_0$. However, in some applications where we want to verify the code by reproducing exact solutions, one must use interpolation. In this chapter, we will use such a problem. diff --git a/chapter2/hyperelasticity.ipynb b/chapter2/hyperelasticity.ipynb index 0cfa9c6b..92633286 100644 --- a/chapter2/hyperelasticity.ipynb +++ b/chapter2/hyperelasticity.ipynb @@ -355,7 +355,7 @@ "id": "nasty-entertainment", "metadata": {}, "source": [ - "Finally, we solve the problem over several time steps, updating the y-component of the traction" + "Finally, we solve the problem over several time steps, updating the z-component of the traction" ] }, { diff --git a/chapter2/hyperelasticity.py b/chapter2/hyperelasticity.py index 2dac5659..2e89cf15 100644 --- a/chapter2/hyperelasticity.py +++ b/chapter2/hyperelasticity.py @@ -180,7 +180,7 @@ def right(x): warped["mag"] = magnitude.x.array # - -# Finally, we solve the problem over several time steps, updating the y-component of the traction +# Finally, we solve the problem over several time steps, updating the z-component of the traction log.set_log_level(log.LogLevel.INFO) tval0 = -1.5 diff --git a/chapter2/linearelasticity.md b/chapter2/linearelasticity.md index ca97184b..f1131b3d 100644 --- a/chapter2/linearelasticity.md +++ b/chapter2/linearelasticity.md @@ -26,11 +26,11 @@ The variational formulation of the PDE consists of forming the inner product of ``` Since $\nabla \cdot \sigma$ contains second-order derivatives of our unknown $u$, we integrate this term by parts ```{math} - -\int_{\Omega}(\nabla \cdot \sigma)\cdot v ~\mathrm{d} x =\int_{\Omega}\sigma : \nabla v ~\mathrm{d}x - \int_{\partial\Omega} (\sigma \cdot n)\cdot v~\mathrm{d}x, + -\int_{\Omega}(\nabla \cdot \sigma)\cdot v ~\mathrm{d} x =\int_{\Omega}\sigma : \nabla v ~\mathrm{d}x - \int_{\partial\Omega} (\sigma \cdot n)\cdot v~\mathrm{d}s, ``` where the colon operator is the inner product between tensors (summed pairwise product of all elements), and $n$ is the outward unit normal at the boundary. The quantity $\sigma \cdot n$ is known as the *traction* or stress vector at the boundary, and often prescribed as a boundary condition. We here assume that it is prescribed on a part $\partial \Omega_T$ of the boundary as $\sigma \cdot n=T$. On the remaining part of the boundary, we assume that the value of the displacement is given as Dirichlet condition (and hence the boundary integral on those boundaries are $0$). We thus obtain ```{math} - \int_{\Omega} \sigma : \nabla v ~\mathrm{d} x = \int_{\Omega} f\cdot v ~\mathrm{d} x + \int_{\partial\Omega_T}Tv~\mathrm{d} s. + \int_{\Omega} \sigma : \nabla v ~\mathrm{d} x = \int_{\Omega} f\cdot v ~\mathrm{d} x + \int_{\partial\Omega_T}T\cdot v~\mathrm{d} s. ``` If we now insert for $\sigma$ its representation with the unknown $u$, we can obtain our variational formulation: Find $u\in V$ such that diff --git a/chapter2/linearelasticity_code.ipynb b/chapter2/linearelasticity_code.ipynb index 07b99e55..d8f54faf 100644 --- a/chapter2/linearelasticity_code.ipynb +++ b/chapter2/linearelasticity_code.ipynb @@ -14,7 +14,7 @@ "- Compute Von Mises stresses\n", "\n", "## Test problem\n", - "As a test example, we will model a clamped beam deformed under its own weigth in 3D. This can be modeled, by setting the right-hand side body force per unit volume to $f=(0,0,-\\rho g)$ with $\\rho$ the density of the beam and $g$ the acceleration of gravity. The beam is box-shaped with length $L$ and has a square cross section of width $W$. we set $u=u_D=(0,0,0)$ at the clamped end, x=0. The rest of the boundary is traction free, that is, we set $T=0$. We start by defining the physical variables used in the program." + "As a test example, we will model a clamped beam deformed under its own weigth in 3D. This can be modeled, by setting the right-hand side body force per unit volume to $f=(0,0,-\\rho g)$ with $\\rho$ the density of the beam and $g$ the acceleration of gravity. The beam is box-shaped with length $L$ and has a square cross section of width $W$. We set $u=u_D=(0,0,0)$ at the clamped end, x=0. The rest of the boundary is traction free, that is, we set $T=0$. We start by defining the physical variables used in the program." ] }, { @@ -160,7 +160,7 @@ "$\\nabla u = \\left(\\frac{\\partial u}{\\partial x}, \\frac{\\partial u}{\\partial y}, \\frac{\\partial u}{\\partial z} \\right)$.\n", "\n", "However, if $u$ is vector valued, the meaning is less clear. Some sources define $\\nabla u$ as a matrix with the elements $\\frac{\\partial u_j}{\\partial x_i}$, while other sources prefer\n", - "$\\frac{\\partial u_i}{\\partial x_j}$. In DOLFINx `grad(u)` is defined as the amtrix with element $\\frac{\\partial u_i}{\\partial x_j}$. However, as it is common in continuum mechanics to use the other definition, `ufl` supplies us with `nabla_grad` for this purpose.\n", + "$\\frac{\\partial u_i}{\\partial x_j}$. In DOLFINx `grad(u)` is defined as the matrix with elements $\\frac{\\partial u_i}{\\partial x_j}$. However, as it is common in continuum mechanics to use the other definition, `ufl` supplies us with `nabla_grad` for this purpose.\n", "```\n", "\n", "## Solve the linear variational problem\n", @@ -1137,7 +1137,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/chapter2/linearelasticity_code.py b/chapter2/linearelasticity_code.py index 5d9679d1..caad8300 100644 --- a/chapter2/linearelasticity_code.py +++ b/chapter2/linearelasticity_code.py @@ -23,7 +23,7 @@ # - Compute Von Mises stresses # # ## Test problem -# As a test example, we will model a clamped beam deformed under its own weigth in 3D. This can be modeled, by setting the right-hand side body force per unit volume to $f=(0,0,-\rho g)$ with $\rho$ the density of the beam and $g$ the acceleration of gravity. The beam is box-shaped with length $L$ and has a square cross section of width $W$. we set $u=u_D=(0,0,0)$ at the clamped end, x=0. The rest of the boundary is traction free, that is, we set $T=0$. We start by defining the physical variables used in the program. +# As a test example, we will model a clamped beam deformed under its own weigth in 3D. This can be modeled, by setting the right-hand side body force per unit volume to $f=(0,0,-\rho g)$ with $\rho$ the density of the beam and $g$ the acceleration of gravity. The beam is box-shaped with length $L$ and has a square cross section of width $W$. We set $u=u_D=(0,0,0)$ at the clamped end, x=0. The rest of the boundary is traction free, that is, we set $T=0$. We start by defining the physical variables used in the program. # Scaled variable import pyvista @@ -99,7 +99,7 @@ def sigma(u): # $\nabla u = \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z} \right)$. # # However, if $u$ is vector valued, the meaning is less clear. Some sources define $\nabla u$ as a matrix with the elements $\frac{\partial u_j}{\partial x_i}$, while other sources prefer -# $\frac{\partial u_i}{\partial x_j}$. In DOLFINx `grad(u)` is defined as the amtrix with element $\frac{\partial u_i}{\partial x_j}$. However, as it is common in continuum mechanics to use the other definition, `ufl` supplies us with `nabla_grad` for this purpose. +# $\frac{\partial u_i}{\partial x_j}$. In DOLFINx `grad(u)` is defined as the matrix with elements $\frac{\partial u_i}{\partial x_j}$. However, as it is common in continuum mechanics to use the other definition, `ufl` supplies us with `nabla_grad` for this purpose. # ``` # # ## Solve the linear variational problem diff --git a/chapter2/navierstokes.md b/chapter2/navierstokes.md index 392cf86e..dd23acc3 100644 --- a/chapter2/navierstokes.md +++ b/chapter2/navierstokes.md @@ -39,7 +39,7 @@ The IPCS scheme involves three steps. First, we compute a *tentative velocity $u -\left\langle \mu \nabla u^{n+\frac{1}{2}}\cdot n, v \right \rangle_{\partial\Omega}= \left\langle f^{n+1}, v \right\rangle. ``` -This notation, suitable for problems wit many terms in the variational formulations, requires some explaination. +This notation, suitable for problems with many terms in the variational formulations, requires some explanation. First, we use the short-hand notation ```{math} \langle v, w \rangle = \int_{\Omega} vw~\mathrm{d}x, \qquad @@ -50,7 +50,7 @@ This allows us to express the variational problem in a more compact way. Second, u^{n+\frac{1}{2}}\approx \frac{u^{n}+ u^{n+1}}{2}. ``` Third, we notice that the variational problem [](ipcs-one) arises from the integration by parts of the term -$langle -\nabla \cdot \sigma, v\rangle$. Just as for the [linear elasticity problem](./linearelasticity.md), we obtain +$\langle -\nabla \cdot \sigma, v\rangle$. Just as for the [linear elasticity problem](./linearelasticity.md), we obtain ```{math} \langle -\nabla \cdot \sigma, v\rangle = \langle \sigma, \epsilon(v) \rangle @@ -67,7 +67,7 @@ which means that we must use the UFL-operator `nabla_grad`. If we use the operat As mentioned in the note in [Linear elasticity implementation](./linearelasticity_code) the usage of `nabla_grad` and `grad` has to be interpreted with care. For the Navier-Stokes equations it is important to consider the term $u\cdot \nabla u$ which should be interpreted as the vector $w$ with elements $w_i=\sum_{j}\left(u_j\frac{\partial}{\partial x_j}\right)u_i = \sum_j u_j\frac{\partial u_i}{\partial x_j}$. This term can be implemented in FEniCSx as either -`grad(u)*u`, since this expression becomes $\sum_j\frac{\partial u_j}{\partial x_j}u_j$, or as `dot(u, nabla_grad(u))` since this +`grad(u)*u`, since this expression becomes $\sum_j\frac{\partial u_i}{\partial x_j}u_j$, or as `dot(u, nabla_grad(u))` since this expression becomes $\sum_i u_i\frac{\partial u_j}{x_i}$. We will use the notation `dot(u, nabla_grad(u))` below since it corresponds more closely to the standard notation $u\cdot \nabla u$. ``` @@ -88,7 +88,7 @@ Taking the divergence and requiring that $\nabla \cdot u^{n+1}=0$ by the Navier- :label: ipcs-tmp - \frac{\rho \nabla\cdot u^*}{\Delta t}+ \nabla^2p^{n+1}-\nabla^2p^n=0, ``` -which is the Poisson problem for the pressure $p^{n+1} resulting in the variational formulation [](ipcs-two). +which is the Poisson problem for the pressure $p^{n+1}$ resulting in the variational formulation [](ipcs-two). Finally, we compute the corrected velocity $u^{n+1}$ from the equation [](ipcs-tmp). Multiplying this equation by a test function $v$, we obtain ```{math} diff --git a/chapter2/ns_code1.ipynb b/chapter2/ns_code1.ipynb index 732f7e71..4b83e79f 100644 --- a/chapter2/ns_code1.ipynb +++ b/chapter2/ns_code1.ipynb @@ -532,7 +532,7 @@ "metadata": {}, "source": [ "## Verification\n", - "As for the previous problems we compute the error at each degree of freedom and the $L^2(\\Omega)$-error. We start with the initial condition $u=(0,0)$. We have not specified the initial condition explicitly, and FEniCSx will initialize all `Function`s including `u_n` and `u_` to zero. Since the exact solution is quadratic, we expect to reach machine precision within finite time. For our implementation, we observe that the error quickly approaches zero, and is of order $10^{-6}$ at $T=10\n", + "As for the previous problems we compute the error at each degree of freedom and the $L^2(\\Omega)$-error. We start with the initial condition $u=(0,0)$. We have not specified the initial condition explicitly, and FEniCSx will initialize all `Function`s including `u_n` and `u_` to zero. Since the exact solution is quadratic, we expect to reach machine precision within finite time. For our implementation, we observe that the error quickly approaches zero, and is of order $10^{-6}$ at $T=10$\n", "\n", "## Visualization of vectors\n", "We have already looked at how to plot higher order functions and vector functions. In this section we will look at how to visualize vector functions with glyphs, instead of warping the mesh." diff --git a/chapter2/ns_code1.py b/chapter2/ns_code1.py index 9a7a2219..a34da692 100644 --- a/chapter2/ns_code1.py +++ b/chapter2/ns_code1.py @@ -316,7 +316,7 @@ def u_exact(x): solver3.destroy() # ## Verification -# As for the previous problems we compute the error at each degree of freedom and the $L^2(\Omega)$-error. We start with the initial condition $u=(0,0)$. We have not specified the initial condition explicitly, and FEniCSx will initialize all `Function`s including `u_n` and `u_` to zero. Since the exact solution is quadratic, we expect to reach machine precision within finite time. For our implementation, we observe that the error quickly approaches zero, and is of order $10^{-6}$ at $T=10 +# As for the previous problems we compute the error at each degree of freedom and the $L^2(\Omega)$-error. We start with the initial condition $u=(0,0)$. We have not specified the initial condition explicitly, and FEniCSx will initialize all `Function`s including `u_n` and `u_` to zero. Since the exact solution is quadratic, we expect to reach machine precision within finite time. For our implementation, we observe that the error quickly approaches zero, and is of order $10^{-6}$ at $T=10$ # # ## Visualization of vectors # We have already looked at how to plot higher order functions and vector functions. In this section we will look at how to visualize vector functions with glyphs, instead of warping the mesh. diff --git a/chapter2/ns_code2.ipynb b/chapter2/ns_code2.ipynb index 95a2bf0b..0d5c1684 100644 --- a/chapter2/ns_code2.ipynb +++ b/chapter2/ns_code2.ipynb @@ -128,7 +128,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by compute the center of mass for each geometrical entitiy." + "To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical entity.\n" ] }, { @@ -557,7 +557,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will also evaluate the pressure at two points, on in front of the obstacle, $(0.15, 0.2)$, and one behind the obstacle, $(0.25, 0.2)$. To do this, we have to find which cell is containing each of the points, so that we can create a linear combination of the local basis functions and coefficients." + "We will also evaluate the pressure at two points, one in front of the obstacle, $(0.15, 0.2)$, and one behind the obstacle, $(0.25, 0.2)$. To do this, we have to find which cell contains each of the points, so that we can create a linear combination of the local basis functions and coefficients.\n" ] }, { @@ -586,7 +586,7 @@ "Other temporal discretization schemes such as the second order backward difference discretization or Crank-Nicholson discretization with Adams-Bashforth linearization are better behaved than our simple backward difference scheme.\n", "```\n", "\n", - "As in the previous example, we create output files for the velocity and pressure and solve the time-dependent problem. As we are solving a time dependent problem with many time steps, we use the `tqdm`-package to visualize the progress. This package can be install with `pip3`." + "As in the previous example, we create output files for the velocity and pressure and solve the time-dependent problem. As we are solving a time dependent problem with many time steps, we use the `tqdm`-package to visualize the progress. This package can be installed with `pip3`.\n" ] }, { @@ -775,7 +775,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/chapter2/ns_code2.py b/chapter2/ns_code2.py index 2fbade61..30c18bf6 100644 --- a/chapter2/ns_code2.py +++ b/chapter2/ns_code2.py @@ -92,7 +92,8 @@ gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid") -# To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by compute the center of mass for each geometrical entitiy. +# To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical entity. +# inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5 inflow, outflow, walls, obstacle = [], [], [], [] @@ -341,7 +342,8 @@ def __call__(self, x): t_u = np.zeros(num_steps, dtype=np.float64) t_p = np.zeros(num_steps, dtype=np.float64) -# We will also evaluate the pressure at two points, on in front of the obstacle, $(0.15, 0.2)$, and one behind the obstacle, $(0.25, 0.2)$. To do this, we have to find which cell is containing each of the points, so that we can create a linear combination of the local basis functions and coefficients. +# We will also evaluate the pressure at two points, one in front of the obstacle, $(0.15, 0.2)$, and one behind the obstacle, $(0.25, 0.2)$. To do this, we have to find which cell contains each of the points, so that we can create a linear combination of the local basis functions and coefficients. +# tree = bb_tree(mesh, mesh.geometry.dim) points = np.array([[0.15, 0.2, 0], [0.25, 0.2, 0]]) @@ -358,7 +360,8 @@ def __call__(self, x): # Other temporal discretization schemes such as the second order backward difference discretization or Crank-Nicholson discretization with Adams-Bashforth linearization are better behaved than our simple backward difference scheme. # ``` # -# As in the previous example, we create output files for the velocity and pressure and solve the time-dependent problem. As we are solving a time dependent problem with many time steps, we use the `tqdm`-package to visualize the progress. This package can be install with `pip3`. +# As in the previous example, we create output files for the velocity and pressure and solve the time-dependent problem. As we are solving a time dependent problem with many time steps, we use the `tqdm`-package to visualize the progress. This package can be installed with `pip3`. +# from pathlib import Path folder = Path("results") diff --git a/chapter3/component_bc.ipynb b/chapter3/component_bc.ipynb index 3e26316d..5ec3dbd1 100644 --- a/chapter3/component_bc.ipynb +++ b/chapter3/component_bc.ipynb @@ -11,6 +11,7 @@ "We will illustrate the problem using a `VectorElement`. However, the method generalizes to any `MixedElement`.\n", "\n", "We will use a slightly modified version of [the linear elasticity demo](./../chapter2/linearelasticity_code), namely\n", + "\n", "$$\n", "-\\nabla \\cdot \\sigma (u) = f\\quad \\text{in } \\Omega,\n", "$$\n", @@ -32,7 +33,7 @@ "$$\n", "We will consider a two dimensional box spanning $[0,L]\\times[0,H]$, where\n", "$\\partial\\Omega_N$ is the left and right side of the beam, $\\partial\\Omega_D$ the bottom of the beam, while $\\partial\\Omega_{Dx}$ is the right side of the beam.\n", - "We will prescribe a displacement $u_x=0$ on the right side of the beam, while the beam is being deformed under its own weight. The sides of the box is traction free." + "We will prescribe a displacement $u_x=0$ on the right side of the beam, while the beam is being deformed under its own weight. The sides of the box are traction free." ] }, { diff --git a/chapter3/component_bc.py b/chapter3/component_bc.py index aa8b8d82..74864625 100644 --- a/chapter3/component_bc.py +++ b/chapter3/component_bc.py @@ -20,6 +20,7 @@ # We will illustrate the problem using a `VectorElement`. However, the method generalizes to any `MixedElement`. # # We will use a slightly modified version of [the linear elasticity demo](./../chapter2/linearelasticity_code), namely +# # $$ # -\nabla \cdot \sigma (u) = f\quad \text{in } \Omega, # $$ @@ -41,7 +42,7 @@ # $$ # We will consider a two dimensional box spanning $[0,L]\times[0,H]$, where # $\partial\Omega_N$ is the left and right side of the beam, $\partial\Omega_D$ the bottom of the beam, while $\partial\Omega_{Dx}$ is the right side of the beam. -# We will prescribe a displacement $u_x=0$ on the right side of the beam, while the beam is being deformed under its own weight. The sides of the box is traction free. +# We will prescribe a displacement $u_x=0$ on the right side of the beam, while the beam is being deformed under its own weight. The sides of the box are traction free. from dolfinx.plot import vtk_mesh import pyvista diff --git a/chapter3/em.ipynb b/chapter3/em.ipynb index 1abe6a46..3f5f1182 100644 --- a/chapter3/em.ipynb +++ b/chapter3/em.ipynb @@ -18,7 +18,7 @@ "We would like to compute the magnetic field $B$ in the iron cylinder, the copper wires, and the surrounding vaccum.\n", "\n", "We start by simplifying the problem to a 2D problem. We can do this by assuming that the cylinder extends far along the z-axis and as a consequence the field is virtually independent of the z-coordinate.\n", - "Next, we consder Maxwell's equation to derive a Poisson equation for the magnetic field (or rather its potential)\n", + "Next, we consider Maxwell's equation to derive a Poisson equation for the magnetic field (or rather its potential)\n", "\n", "$$\n", "\\nabla \\cdot D = \\rho,\n", @@ -58,14 +58,14 @@ "\\lim_{\\vert(x,y)\\vert\\to \\infty}A_z = 0.\n", "$$\n", "\n", - "Since we cannot solve the problem on an infinite domain, we will truncate the domain using a large disk, and set $A_z=0$ on the boundary. The current $J_z$ is set to $+1$A in the interior set of the circles (copper-wire cross sections) and to $-1$ A in the exteriror set of circles in the cross section figure.\n", + "Since we cannot solve the problem on an infinite domain, we will truncate the domain using a large disk, and set $A_z=0$ on the boundary. The current $J_z$ is set to $+1$A in the interior set of the circles (copper-wire cross sections) and to $-1$ A in the exterior set of circles in the cross section figure.\n", "Once the magnetic field vector potential has been computed, we can compute the magnetic field $B=B(x,y)$ by\n", "\n", "$$\n", " B(x,y)=\\left(\\frac{\\partial A_z}{\\partial y}, - \\frac{\\partial A_z}{\\partial x} \\right).\n", "$$\n", "\n", - "The weak formulation is easily obtained by multiplication of a test function $v$, followed by integration by parts, where all boundary integrals vanishes due to the Dirichlet condition, we obtain $a(A_z,v)=L(v)$ with\n", + "The weak formulation is easily obtained by multiplication of a test function $v$, followed by integration by parts, where all boundary integrals vanish due to the Dirichlet condition, we obtain $a(A_z,v)=L(v)$ with\n", "\n", "$$\n", "a(A_z, v)=\\int_\\Omega \\mu^{-1}\\nabla A_z \\cdot \\nabla v ~\\mathrm{d}x,\n", @@ -1098,7 +1098,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, we define the discontinous functions for the permability $\\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains)" + "Next, we define the discontinous functions for the permeability $\\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains)\n" ] }, { @@ -1480,7 +1480,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/chapter3/em.py b/chapter3/em.py index 23926d12..058b615b 100644 --- a/chapter3/em.py +++ b/chapter3/em.py @@ -27,7 +27,7 @@ # We would like to compute the magnetic field $B$ in the iron cylinder, the copper wires, and the surrounding vaccum. # # We start by simplifying the problem to a 2D problem. We can do this by assuming that the cylinder extends far along the z-axis and as a consequence the field is virtually independent of the z-coordinate. -# Next, we consder Maxwell's equation to derive a Poisson equation for the magnetic field (or rather its potential) +# Next, we consider Maxwell's equation to derive a Poisson equation for the magnetic field (or rather its potential) # # $$ # \nabla \cdot D = \rho, @@ -67,14 +67,14 @@ # \lim_{\vert(x,y)\vert\to \infty}A_z = 0. # $$ # -# Since we cannot solve the problem on an infinite domain, we will truncate the domain using a large disk, and set $A_z=0$ on the boundary. The current $J_z$ is set to $+1$A in the interior set of the circles (copper-wire cross sections) and to $-1$ A in the exteriror set of circles in the cross section figure. +# Since we cannot solve the problem on an infinite domain, we will truncate the domain using a large disk, and set $A_z=0$ on the boundary. The current $J_z$ is set to $+1$A in the interior set of the circles (copper-wire cross sections) and to $-1$ A in the exterior set of circles in the cross section figure. # Once the magnetic field vector potential has been computed, we can compute the magnetic field $B=B(x,y)$ by # # $$ # B(x,y)=\left(\frac{\partial A_z}{\partial y}, - \frac{\partial A_z}{\partial x} \right). # $$ # -# The weak formulation is easily obtained by multiplication of a test function $v$, followed by integration by parts, where all boundary integrals vanishes due to the Dirichlet condition, we obtain $a(A_z,v)=L(v)$ with +# The weak formulation is easily obtained by multiplication of a test function $v$, followed by integration by parts, where all boundary integrals vanish due to the Dirichlet condition, we obtain $a(A_z,v)=L(v)$ with # # $$ # a(A_z, v)=\int_\Omega \mu^{-1}\nabla A_z \cdot \nabla v ~\mathrm{d}x, @@ -225,7 +225,8 @@ cell_tag_fig = plotter.screenshot("cell_tags.png") -# Next, we define the discontinous functions for the permability $\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains) +# Next, we define the discontinous functions for the permeability $\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains) +# # + diff --git a/chapter3/multiple_dirichlet.ipynb b/chapter3/multiple_dirichlet.ipynb index c305aa65..5caa84c5 100644 --- a/chapter3/multiple_dirichlet.ipynb +++ b/chapter3/multiple_dirichlet.ipynb @@ -6,7 +6,7 @@ "source": [ "# Setting multiple Dirichlet condition\n", "\n", - "In the previous section, we used a single function to $u_d$ to setting Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functions, one for each subdomain of the boundary. We consider a similar example to [the previous example](./neumann_dirichlet_code) and redefine it consist of two Dirichlet boundary conditions\n", + "In the previous section, we used a single function for $u_D$ to set Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functions, one for each subdomain of the boundary. We consider a similar example to [the previous example](./neumann_dirichlet_code) and redefine it to consist of two Dirichlet boundary conditions\n", "\n", "$$\n", "-\\nabla^2 u =f \\quad \\text{in } \\Omega,\n", diff --git a/chapter3/multiple_dirichlet.py b/chapter3/multiple_dirichlet.py index 2d0e243f..6b8da4ed 100644 --- a/chapter3/multiple_dirichlet.py +++ b/chapter3/multiple_dirichlet.py @@ -15,7 +15,7 @@ # # Setting multiple Dirichlet condition # -# In the previous section, we used a single function to $u_d$ to setting Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functions, one for each subdomain of the boundary. We consider a similar example to [the previous example](./neumann_dirichlet_code) and redefine it consist of two Dirichlet boundary conditions +# In the previous section, we used a single function for $u_D$ to set Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functions, one for each subdomain of the boundary. We consider a similar example to [the previous example](./neumann_dirichlet_code) and redefine it to consist of two Dirichlet boundary conditions # # $$ # -\nabla^2 u =f \quad \text{in } \Omega, diff --git a/chapter3/neumann_dirichlet_code.ipynb b/chapter3/neumann_dirichlet_code.ipynb index f043b527..3f31e727 100644 --- a/chapter3/neumann_dirichlet_code.ipynb +++ b/chapter3/neumann_dirichlet_code.ipynb @@ -115,7 +115,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fullfils this condition." + "Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fulfill this condition." ] }, { diff --git a/chapter3/neumann_dirichlet_code.py b/chapter3/neumann_dirichlet_code.py index f4adc919..a29b6683 100644 --- a/chapter3/neumann_dirichlet_code.py +++ b/chapter3/neumann_dirichlet_code.py @@ -114,7 +114,7 @@ # - -# Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fullfils this condition. +# Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fulfill this condition. # + def u_exact(x): diff --git a/chapter3/robin_neumann_dirichlet.ipynb b/chapter3/robin_neumann_dirichlet.ipynb index 27307b78..1b9d9dd7 100644 --- a/chapter3/robin_neumann_dirichlet.ipynb +++ b/chapter3/robin_neumann_dirichlet.ipynb @@ -15,7 +15,7 @@ "- $\\Gamma_N$ for Neumann conditions: $-\\kappa \\frac{\\partial u}{\\partial n}=g_j \\text{ on } \\Gamma_N^j$ where $\\Gamma_N=\\Gamma_N^0\\cup \\Gamma_N^1 \\cup \\dots$.\n", "- $\\Gamma_R$ for Robin conditions: $-\\kappa \\frac{\\partial u}{\\partial n}=r(u-s)$\n", "\n", - "where $r$ and $s$ are specified functions. The Robin condition is most often used to model heat transfer to the surroundings and arise naturally from Newton's cooling law.\n", + "where $r$ and $s$ are specified functions. The Robin condition is most often used to model heat transfer to the surroundings and arises naturally from Newton's cooling law.\n", "In that case, $r$ is a heat transfer coefficient, and $s$ is the temperature of the surroundings. \n", "Both can be space and time-dependent. The Robin conditions apply at some parts $\\Gamma_R^0,\\Gamma_R^1,\\dots$, of the boundary:\n", "\n", @@ -125,12 +125,12 @@ " u_D=u_{ex}=1+x^2+2y^2\n", "$$\n", "$$\n", - " g_0=\\left.\\frac{\\partial u_{ex}}{y}\\right\\vert_{y=1}=4y\\vert_{y=1}=-4\n", + " g_0=-\\left.\\frac{\\partial u_{ex}}{\\partial y}\\right\\vert_{y=1}=-4y\\vert_{y=1}=-4\n", "$$\n", "\n", "The Robin condition can be specified in many ways. As\n", - "$-\\left.\\frac{\\partial u_{ex}}{n}\\right\\vert_{x=0}=\\left.\\frac{\\partial u_{ex}}{\\partial x}\\right\\vert_{x=0}=2x=0,$\n", - "we can specify $r\\neq 0$ arbitrarly and $s=u_{ex}$. We choose $r=1000$.\n", + "$-\\left.\\frac{\\partial u_{ex}}{\\partial n}\\right\\vert_{y=0}=\\left.\\frac{\\partial u_{ex}}{\\partial y}\\right\\vert_{y=0}=4y=0,$\n", + "we can specify $r\\neq 0$ arbitrarily and $s=u_{ex}$. We choose $r=1000$.\n", "We can now create all the necessary variable definitions and the traditional part of the variational form." ] }, @@ -352,7 +352,7 @@ "metadata": {}, "source": [ "## Verification\n", - "As for the previous problems, we compute the error of our computed solution and compare it to the analytical solution." + "As for the previous problems, we compute the error of our computed solution and compare it with the analytical solution." ] }, { diff --git a/chapter3/robin_neumann_dirichlet.py b/chapter3/robin_neumann_dirichlet.py index 6a0e9787..1668e6fb 100644 --- a/chapter3/robin_neumann_dirichlet.py +++ b/chapter3/robin_neumann_dirichlet.py @@ -24,7 +24,7 @@ # - $\Gamma_N$ for Neumann conditions: $-\kappa \frac{\partial u}{\partial n}=g_j \text{ on } \Gamma_N^j$ where $\Gamma_N=\Gamma_N^0\cup \Gamma_N^1 \cup \dots$. # - $\Gamma_R$ for Robin conditions: $-\kappa \frac{\partial u}{\partial n}=r(u-s)$ # -# where $r$ and $s$ are specified functions. The Robin condition is most often used to model heat transfer to the surroundings and arise naturally from Newton's cooling law. +# where $r$ and $s$ are specified functions. The Robin condition is most often used to model heat transfer to the surroundings and arises naturally from Newton's cooling law. # In that case, $r$ is a heat transfer coefficient, and $s$ is the temperature of the surroundings. # Both can be space and time-dependent. The Robin conditions apply at some parts $\Gamma_R^0,\Gamma_R^1,\dots$, of the boundary: # @@ -119,12 +119,12 @@ # u_D=u_{ex}=1+x^2+2y^2 # $$ # $$ -# g_0=\left.\frac{\partial u_{ex}}{y}\right\vert_{y=1}=4y\vert_{y=1}=-4 +# g_0=-\left.\frac{\partial u_{ex}}{\partial y}\right\vert_{y=1}=-4y\vert_{y=1}=-4 # $$ # # The Robin condition can be specified in many ways. As -# $-\left.\frac{\partial u_{ex}}{n}\right\vert_{x=0}=\left.\frac{\partial u_{ex}}{\partial x}\right\vert_{x=0}=2x=0,$ -# we can specify $r\neq 0$ arbitrarly and $s=u_{ex}$. We choose $r=1000$. +# $-\left.\frac{\partial u_{ex}}{\partial n}\right\vert_{y=0}=\left.\frac{\partial u_{ex}}{\partial y}\right\vert_{y=0}=4y=0,$ +# we can specify $r\neq 0$ arbitrarily and $s=u_{ex}$. We choose $r=1000$. # We can now create all the necessary variable definitions and the traditional part of the variational form. u_ex = lambda x: 1 + x[0]**2 + 2*x[1]**2 @@ -244,7 +244,7 @@ def type(self): # - # ## Verification -# As for the previous problems, we compute the error of our computed solution and compare it to the analytical solution. +# As for the previous problems, we compute the error of our computed solution and compare it with the analytical solution. # + # Compute L2 error and error at nodes diff --git a/chapter3/subdomains.ipynb b/chapter3/subdomains.ipynb index c1e42bcc..99218d5b 100644 --- a/chapter3/subdomains.ipynb +++ b/chapter3/subdomains.ipynb @@ -7,7 +7,7 @@ "# Defining subdomains for different materials\n", "Author: Jørgen S. Dokken\n", "\n", - "Solving PDEs in domains made up of different materials is frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function.\n", + "Solving PDEs in domains made up of different materials is a frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function.\n", "Such a function can be created over any mesh in the following way\n", "## Subdomains on built-in meshes" ] @@ -68,7 +68,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that both fucntion uses a $\\leq$ or $\\geq$, as FEniCSx will evaluate each cell at all of the vertices, and thus for has to return `True` for all vertices align with the interface to be marked properly.\n", + "Note that both functions use a $\\leq$ or $\\geq$, as FEniCSx will evaluate each cell at all of the vertices, and thus has to return `True` for all vertices aligned with the interface to be marked properly.\n", "\n", "We will solve a variable-coefficient extension of the Poisson equation\n", "\n", @@ -100,7 +100,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the previous code block, we found which cells (triangular elements) which satisfies the condition for being in $\\Omega_0, \\Omega_1$. As the $DG-0$ function contain only one degree of freedom per mesh, there is a one to one mapping between the cell indicies and the degrees of freedom. We let $\\kappa=\\begin{cases}\n", + "In the previous code block, we found which cells (triangular elements) satisfy the condition for being in $\\Omega_0, \\Omega_1$. As the $DG-0$ function contains only one degree of freedom per cell, there is a one to one mapping between the cell indicies and the degrees of freedom. We let $\\kappa=\\begin{cases}\n", "1 &\\text{if } x\\in\\Omega_0\\\\\n", "0.1& \\text{if } x\\in\\Omega_1\\\\\n", "\\end{cases}$" @@ -220,7 +220,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We clearly observe different behavior in the two regions, whose both has the same Dirichlet boundary condition on the left side, where $x=0$." + "We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$." ] }, { diff --git a/chapter3/subdomains.py b/chapter3/subdomains.py index ba9bbd6c..4f6326a3 100644 --- a/chapter3/subdomains.py +++ b/chapter3/subdomains.py @@ -16,7 +16,7 @@ # # Defining subdomains for different materials # Author: Jørgen S. Dokken # -# Solving PDEs in domains made up of different materials is frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function. +# Solving PDEs in domains made up of different materials is a frequently encountered task. In FEniCSx, we handle these problems by defining a Discontinous cell-wise constant function. # Such a function can be created over any mesh in the following way # ## Subdomains on built-in meshes @@ -61,7 +61,7 @@ def Omega_1(x): # - -# Note that both fucntion uses a $\leq$ or $\geq$, as FEniCSx will evaluate each cell at all of the vertices, and thus for has to return `True` for all vertices align with the interface to be marked properly. +# Note that both functions use a $\leq$ or $\geq$, as FEniCSx will evaluate each cell at all of the vertices, and thus has to return `True` for all vertices aligned with the interface to be marked properly. # # We will solve a variable-coefficient extension of the Poisson equation # @@ -81,7 +81,7 @@ def Omega_1(x): cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0) cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1) -# In the previous code block, we found which cells (triangular elements) which satisfies the condition for being in $\Omega_0, \Omega_1$. As the $DG-0$ function contain only one degree of freedom per mesh, there is a one to one mapping between the cell indicies and the degrees of freedom. We let $\kappa=\begin{cases} +# In the previous code block, we found which cells (triangular elements) satisfy the condition for being in $\Omega_0, \Omega_1$. As the $DG-0$ function contains only one degree of freedom per cell, there is a one to one mapping between the cell indicies and the degrees of freedom. We let $\kappa=\begin{cases} # 1 &\text{if } x\in\Omega_0\\ # 0.1& \text{if } x\in\Omega_1\\ # \end{cases}$ @@ -135,7 +135,7 @@ def Omega_1(x): figure = p2.screenshot("subdomains_structured2.png") -# We clearly observe different behavior in the two regions, whose both has the same Dirichlet boundary condition on the left side, where $x=0$. +# We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$. # ## Interpolation with Python-function # As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation: diff --git a/chapter4/compiler_parameters.ipynb b/chapter4/compiler_parameters.ipynb index 437f07bc..d7a8c617 100644 --- a/chapter4/compiler_parameters.ipynb +++ b/chapter4/compiler_parameters.ipynb @@ -11,12 +11,12 @@ "In this chapter, we will explore how to optimize and inspect the integration kernels used in DOLFINx.\n", "As we have seen in the previous demos, DOLFINx uses the [Unified form language](https://github.com/FEniCS/ufl/) to describe variational problems.\n", "\n", - "These descriptions has to be translated in to code for assembling the right and left hand side of the discrete variational problem.\n", + "These descriptions have to be translated into code for assembling the right and left hand side of the discrete variational problem.\n", "\n", "DOLFINx uses [ffcx](https://github.com/FEniCS/ffcx/) to generate efficient C code assembling the element matrices.\n", - "This C code is in turned compiled using [CFFI](https://cffi.readthedocs.io/en/latest/), and we can specify a variety of compile options.\n", + "This C code is in turn compiled using [CFFI](https://cffi.readthedocs.io/en/latest/), and we can specify a variety of compile options.\n", "\n", - "We start by specifying the current directory as the place to place the generated C files, we obtain the current directory using pathlib" + "We start by specifying the current directory as the location to place the generated C files, we obtain the current directory using pathlib" ] }, { @@ -89,7 +89,7 @@ "id": "93afb3b1", "metadata": {}, "source": [ - "We start by considering the different levels of optimization the C compiled can use on the optimized code. A list of optimization options and explainations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)" + "We start by considering the different levels of optimization that the C compiler can use on the optimized code. A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)" ] }, { diff --git a/chapter4/compiler_parameters.py b/chapter4/compiler_parameters.py index 16e676ba..756b11dd 100644 --- a/chapter4/compiler_parameters.py +++ b/chapter4/compiler_parameters.py @@ -19,12 +19,12 @@ # In this chapter, we will explore how to optimize and inspect the integration kernels used in DOLFINx. # As we have seen in the previous demos, DOLFINx uses the [Unified form language](https://github.com/FEniCS/ufl/) to describe variational problems. # -# These descriptions has to be translated in to code for assembling the right and left hand side of the discrete variational problem. +# These descriptions have to be translated into code for assembling the right and left hand side of the discrete variational problem. # # DOLFINx uses [ffcx](https://github.com/FEniCS/ffcx/) to generate efficient C code assembling the element matrices. -# This C code is in turned compiled using [CFFI](https://cffi.readthedocs.io/en/latest/), and we can specify a variety of compile options. +# This C code is in turn compiled using [CFFI](https://cffi.readthedocs.io/en/latest/), and we can specify a variety of compile options. # -# We start by specifying the current directory as the place to place the generated C files, we obtain the current directory using pathlib +# We start by specifying the current directory as the location to place the generated C files, we obtain the current directory using pathlib # + import matplotlib.pyplot as plt @@ -67,7 +67,7 @@ def compile_form(space: str, degree: int, jit_options: Dict): # - -# We start by considering the different levels of optimization the C compiled can use on the optimized code. A list of optimization options and explainations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) +# We start by considering the different levels of optimization that the C compiler can use on the optimized code. A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) optimization_options = ["-O1", "-O2", "-O3", "-Ofast"] diff --git a/chapter4/convergence.ipynb b/chapter4/convergence.ipynb index 1b83ac5d..bfc532ba 100644 --- a/chapter4/convergence.ipynb +++ b/chapter4/convergence.ipynb @@ -70,7 +70,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now, we can compute the error between the analyical solution `u_ex=u_ufl(x)` and approximated solution `uh`. A natural choice might seem to compute `(u_ex-uh)**2*ufl.dx`." + "Now, we can compute the error between the analyical solution `u_ex=u_ufl(x)` and the approximated solution `uh`. A natural choice might seem to compute `(u_ex-uh)**2*ufl.dx`." ] }, { @@ -99,7 +99,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Sometimes it is of interest to compute the error fo the gradient field, $\\vert\\vert \\nabla(u_e-u_h)\\vert\\vert$, often referred to as the $H_0^1$-nrom of the error, this can be expressed as" + "Sometimes it is of interest to compute the error fo the gradient field, $\\vert\\vert \\nabla(u_e-u_h)\\vert\\vert$, often referred to as the $H_0^1$-norm of the error, this can be expressed as" ] }, { @@ -214,11 +214,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "If we assume that $E_i$ is of the form $E_i=Ch_i^r$, with unknown constants $C$ and $r$, we can compare two consecqutive experiments, $E_{i-1}= Ch_{i-1}^r$ and $E_i=Ch_i^r$, and solve for $r$:\n", + "If we assume that $E_i$ is of the form $E_i=Ch_i^r$, with unknown constants $C$ and $r$, we can compare two consecutive experiments, $E_{i-1}= Ch_{i-1}^r$ and $E_i=Ch_i^r$, and solve for $r$:\n", "```{math}\n", "r=\\frac{\\ln(E_i/E_{i-1})}{\\ln(h_i/h_{i-1})}\n", "```\n", - "The $r$ values should approac the expected convergence rate (which is typically the polynomial degree + 1 for the $L^2$-error.) as $i$ increases. This can be written compactly using `numpy`." + "The $r$ values should approach the expected convergence rate (which is typically the polynomial degree + 1 for the $L^2$-error.) as $i$ increases. This can be written compactly using `numpy`." ] }, { diff --git a/chapter4/convergence.py b/chapter4/convergence.py index e39c0dc2..e32fe655 100644 --- a/chapter4/convergence.py +++ b/chapter4/convergence.py @@ -67,7 +67,7 @@ def solve_poisson(N=10, degree=1): # - -# Now, we can compute the error between the analyical solution `u_ex=u_ufl(x)` and approximated solution `uh`. A natural choice might seem to compute `(u_ex-uh)**2*ufl.dx`. +# Now, we can compute the error between the analyical solution `u_ex=u_ufl(x)` and the approximated solution `uh`. A natural choice might seem to compute `(u_ex-uh)**2*ufl.dx`. uh, u_ex = solve_poisson(10) comm = uh.function_space.mesh.comm @@ -76,7 +76,7 @@ def solve_poisson(N=10, degree=1): if comm.rank == 0: print(f"L2-error: {E:.2e}") -# Sometimes it is of interest to compute the error fo the gradient field, $\vert\vert \nabla(u_e-u_h)\vert\vert$, often referred to as the $H_0^1$-nrom of the error, this can be expressed as +# Sometimes it is of interest to compute the error fo the gradient field, $\vert\vert \nabla(u_e-u_h)\vert\vert$, often referred to as the $H_0^1$-norm of the error, this can be expressed as eh = uh - u_ex error_H10 = form(dot(grad(eh), grad(eh)) * dx) @@ -137,11 +137,11 @@ def error_L2(uh, u_ex, degree_raise=3): if comm.rank == 0: print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}") -# If we assume that $E_i$ is of the form $E_i=Ch_i^r$, with unknown constants $C$ and $r$, we can compare two consecqutive experiments, $E_{i-1}= Ch_{i-1}^r$ and $E_i=Ch_i^r$, and solve for $r$: +# If we assume that $E_i$ is of the form $E_i=Ch_i^r$, with unknown constants $C$ and $r$, we can compare two consecutive experiments, $E_{i-1}= Ch_{i-1}^r$ and $E_i=Ch_i^r$, and solve for $r$: # ```{math} # r=\frac{\ln(E_i/E_{i-1})}{\ln(h_i/h_{i-1})} # ``` -# The $r$ values should approac the expected convergence rate (which is typically the polynomial degree + 1 for the $L^2$-error.) as $i$ increases. This can be written compactly using `numpy`. +# The $r$ values should approach the expected convergence rate (which is typically the polynomial degree + 1 for the $L^2$-error.) as $i$ increases. This can be written compactly using `numpy`. rates = np.log(Es[1:] / Es[:-1]) / np.log(hs[1:] / hs[:-1]) if comm.rank == 0: diff --git a/chapter4/newton-solver.ipynb b/chapter4/newton-solver.ipynb index da40c87c..83ceb616 100644 --- a/chapter4/newton-solver.ipynb +++ b/chapter4/newton-solver.ipynb @@ -305,7 +305,7 @@ "metadata": {}, "source": [ "## Visualization of Newton iterations\n", - "We next look at the evolution of the solutions and the error of the solution when compared to the two exact roots of the problem." + "We next look at the evolution of the solution and the error of the solution when compared to the two exact roots of the problem." ] }, { diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index b0400815..af4b12a9 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -148,7 +148,7 @@ def root_1(x): print(f"Final residual {L.norm(0)}") # ## Visualization of Newton iterations -# We next look at the evolution of the solutions and the error of the solution when compared to the two exact roots of the problem. +# We next look at the evolution of the solution and the error of the solution when compared to the two exact roots of the problem. # + # Plot solution for each of the iterations diff --git a/chapter4/solvers.ipynb b/chapter4/solvers.ipynb index f09d9132..7d095c01 100644 --- a/chapter4/solvers.ipynb +++ b/chapter4/solvers.ipynb @@ -7,7 +7,7 @@ "# Solver configuration\n", "Author: Jørgen S. Dokken\n", "\n", - "In this section, we will go through how to specify what linear algebra solver we would like to use to solve our PDEs, as well as how to verify the implemenation by considering convergence rates.\n", + "In this section, we will go through how to specify what linear algebra solver we would like to use to solve our PDEs, as well as how to verify the implementation by considering convergence rates.\n", "\n", "$$\n", "-\\Delta u = f \\text{ in } \\Omega\n", @@ -205,9 +205,9 @@ "source": [ "This is a very robust and simple method, and is the recommended method up to a few thousand unknowns and can be efficiently used for many 2D and smaller 3D problems. However, sparse LU decomposition quickly becomes slow, as for a $N\\times N$-matrix the number of floating point operations scales as $\\sim (2/3)N^3$.\n", "\n", - "For large problems, we instead need to use an iterative method which are faster and require less memory.\n", + "For large problems, we instead need to use an iterative method which is faster and requires less memory.\n", "## Choosing a linear solver and preconditioner\n", - "As the Poisson equation results in a symmetric, positive definite system matrix, the optimal Krylov solver is the conjugate gradient (Lagrange) method. The default preconditioner is the incomplete LU factorization (ILU), which is a popular and robous overall preconditioner. We can change the preconditioner by setting `\"pc_type\"` to some of the other preconditioners in petsc, which you can find in at [PETSc KSP solvers](https://petsc.org/release/manual/ksp/#tab-kspdefaults) and [PETSc preconditioners](https://petsc.org/release/manual/ksp/#tab-pcdefaults).\n", + "As the Poisson equation results in a symmetric, positive definite system matrix, the optimal Krylov solver is the conjugate gradient (Lagrange) method. The default preconditioner is the incomplete LU factorization (ILU), which is a popular and robust overall preconditioner. We can change the preconditioner by setting `\"pc_type\"` to some of the other preconditioners in petsc, which you can find at [PETSc KSP solvers](https://petsc.org/release/manual/ksp/#tab-kspdefaults) and [PETSc preconditioners](https://petsc.org/release/manual/ksp/#tab-pcdefaults).\n", "You can set any option in `PETSc` through the `petsc_options` input, such as the absolute tolerance (`\"ksp_atol\"`), relative tolerance (`\"ksp_rtol\"`) and maximum number of iterations (`\"ksp_max_it\"`)." ] }, @@ -293,7 +293,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For non-symmetrix problems, a Krylov solver for non-symmetrix systems, such as GMRES is a better." + "For non-symmetric problems, a Krylov solver for non-symmetric systems, such as GMRES is better." ] }, { diff --git a/chapter4/solvers.py b/chapter4/solvers.py index 6194d2c2..7866099c 100644 --- a/chapter4/solvers.py +++ b/chapter4/solvers.py @@ -16,7 +16,7 @@ # # Solver configuration # Author: Jørgen S. Dokken # -# In this section, we will go through how to specify what linear algebra solver we would like to use to solve our PDEs, as well as how to verify the implemenation by considering convergence rates. +# In this section, we will go through how to specify what linear algebra solver we would like to use to solve our PDEs, as well as how to verify the implementation by considering convergence rates. # # $$ # -\Delta u = f \text{ in } \Omega @@ -87,9 +87,9 @@ def u_ex(mod): # This is a very robust and simple method, and is the recommended method up to a few thousand unknowns and can be efficiently used for many 2D and smaller 3D problems. However, sparse LU decomposition quickly becomes slow, as for a $N\times N$-matrix the number of floating point operations scales as $\sim (2/3)N^3$. # -# For large problems, we instead need to use an iterative method which are faster and require less memory. +# For large problems, we instead need to use an iterative method which is faster and requires less memory. # ## Choosing a linear solver and preconditioner -# As the Poisson equation results in a symmetric, positive definite system matrix, the optimal Krylov solver is the conjugate gradient (Lagrange) method. The default preconditioner is the incomplete LU factorization (ILU), which is a popular and robous overall preconditioner. We can change the preconditioner by setting `"pc_type"` to some of the other preconditioners in petsc, which you can find in at [PETSc KSP solvers](https://petsc.org/release/manual/ksp/#tab-kspdefaults) and [PETSc preconditioners](https://petsc.org/release/manual/ksp/#tab-pcdefaults). +# As the Poisson equation results in a symmetric, positive definite system matrix, the optimal Krylov solver is the conjugate gradient (Lagrange) method. The default preconditioner is the incomplete LU factorization (ILU), which is a popular and robust overall preconditioner. We can change the preconditioner by setting `"pc_type"` to some of the other preconditioners in petsc, which you can find at [PETSc KSP solvers](https://petsc.org/release/manual/ksp/#tab-kspdefaults) and [PETSc preconditioners](https://petsc.org/release/manual/ksp/#tab-pcdefaults). # You can set any option in `PETSc` through the `petsc_options` input, such as the absolute tolerance (`"ksp_atol"`), relative tolerance (`"ksp_rtol"`) and maximum number of iterations (`"ksp_max_it"`). cg_problem = LinearProblem(a, L, bcs=bcs, @@ -102,7 +102,7 @@ def u_ex(mod): for line in solver_output.readlines(): print(line) -# For non-symmetrix problems, a Krylov solver for non-symmetrix systems, such as GMRES is a better. +# For non-symmetric problems, a Krylov solver for non-symmetric systems, such as GMRES is better. gmres_problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "gmres", "ksp_rtol": 1e-6, "ksp_atol": 1e-10, "ksp_max_it": 1000, "pc_type": "none"})