diff --git a/FEMII/00_intro_to_fenics.ipynb b/FEMII/00_intro_to_fenics.ipynb
index 70ac73c..68a6156 100644
--- a/FEMII/00_intro_to_fenics.ipynb
+++ b/FEMII/00_intro_to_fenics.ipynb
@@ -1,699 +1,1372 @@
{
"cells": [
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# IPython magic to import matplotlib and plot inline\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"# Introduction to FEniCS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[FEniCS](https://fenicsproject.org/) is a open-source computing platform for solving partial differential equations (PDEs). \n",
"The FEniCS Project consists of a number of components with DOLFIN and UFL providing the main user interface. For Detailed documentation of the FEniCS programming interface and its components, see the links on https://fenics.readthedocs.io/en/latest/.\n",
"\n",
"\n",
"[The FEniCS Book](https://link.springer.com/book/10.1007%2F978-3-642-23099-8) provides an overview of the finite element method, its implementation in FEniCS and application examples.\n",
"[The FEniCS Tutorial](https://fenicsproject.org/tutorial/) is a collection of example problems for learning FEniCS 'by examples'. \n",
"Source code for the examples used this book can be found [here](https://github.com/hplgit/fenics-tutorial/tree/master/src/vol1/python).\n",
"Other examples can be found in the [DOLFIN documentation](https://fenicsproject.org/docs/dolfin/latest/python/demos.html).\n",
+ "A recent FEniCS hands-on tutorial is available [here](https://fenics-handson.readthedocs.io/en/latest/index.html).\n",
"\n",
- "Be aware that FEniCS is under development and unfortunately some class and function names of the current version may differ from those in the tutorials...\n",
+ "Be aware that some class and function names of the current version of FEniCS may differ from those in the tutorials...\n",
"\n",
"This notebook assumes FEniCS/Dolfin version `2019.1.0`. \n",
"\n",
"You are running this notebook with version:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Your Dolfin version: 2019.1.0\n"
+ ]
+ }
+ ],
"source": [
"from dolfin import __version__ as dolfin_version\n",
"print(\"Your Dolfin version: \", dolfin_version)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting Started"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Follow the instructions in the *README.md* file in this repository to set up a running fenics installation.\n",
"- Skim read section 1 of [The FEniCS Tutorial](https://fenicsproject.org/pub/tutorial/html/._ftut1003.html#ch:prelim)\n",
"- Work through this notebook. This will give you a quick overview of the basic building blocks of any FEniCS program, and introduce an example how PDE problems can be solved in FEniCS. \n",
"- Then, read and follow the examples of section 2 in [The FEniCS Tutorial](https://fenicsproject.org/pub/tutorial/html/._ftut1004.html#ch:fundamentals). You will already be familiar with many of the concepts, but they will be explained in more detail.\n",
"\n",
"After this, you are ready to explore any of the other examples in the tutorial.\n",
"Particularly relevant for your projects are:\n",
"- [Linear Elasticity](https://fenicsproject.org/pub/tutorial/html/._ftut1008.html#ftut:elast) example in the tutorial.\n",
- "- [Hyperelasticity](https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/hyperelasticity/python/documentation.html) from the Dolfin reference. Note that this code is based on an outdated version of Dolfing/FEniCS, so a few adapatations will be necessary to make it run.\n",
+ "- [Hyperelasticity](https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/hyperelasticity/python/documentation.html) from the Dolfin reference. Note that this code is based on an outdated version of Dolfing/FEniCS, so a few adapatations will be necessary to make it run. A more recent (also more complex) example is available [here](https://fenics-handson.readthedocs.io/en/latest/elasticity/doc.html).\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importing FEniCS libraries"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from fenics import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mesh Creation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"FEniCS provides various [built-in meshes](https://fenicsproject.org/docs/dolfin/latest/python/demos/built-in-meshes/demo_built-in-meshes.py.html) for simple standard geometries.\n",
"Meshes for more complex geometries can be generated via [fenics-mshr](https://fenics.readthedocs.io/projects/mshr/en/latest/) / [CGAL](https://www.cgal.org)."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 4,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ " \n",
+ " \n",
+ " "
+ ],
+ "text/html": [
+ "\n",
+ "\n",
+ " \n",
+ " FEniCS/DOLFIN X3DOM plot \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"mesh_2d = RectangleMesh(Point(0,0),Point(5,5), nx=3, ny=3)\n",
"# plotting via ipython 3d rendering \n",
"mesh_2d\n",
"# or via matplotlib in 2d\n",
"#plot(mesh_2d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "A mesh consists of various 'entities' as detailed [here](https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/cpp/mesh/Mesh.html)."
+ "A mesh consists of various 'entities' as detailed [here](https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/cpp/mesh/Mesh.html)."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Number of cells: 18\n",
+ "Number of faces: 18\n",
+ "Number of edges: 33\n",
+ "Number of vertices: 16\n",
+ "Coordinates: \n",
+ " [[0. 0. ]\n",
+ " [1.66666667 0. ]\n",
+ " [3.33333333 0. ]\n",
+ " [5. 0. ]\n",
+ " [0. 1.66666667]\n",
+ " [1.66666667 1.66666667]\n",
+ " [3.33333333 1.66666667]\n",
+ " [5. 1.66666667]\n",
+ " [0. 3.33333333]\n",
+ " [1.66666667 3.33333333]\n",
+ " [3.33333333 3.33333333]\n",
+ " [5. 3.33333333]\n",
+ " [0. 5. ]\n",
+ " [1.66666667 5. ]\n",
+ " [3.33333333 5. ]\n",
+ " [5. 5. ]]\n"
+ ]
+ }
+ ],
"source": [
"print(\"Number of cells: \", mesh_2d.num_cells())\n",
"print(\"Number of faces: \", mesh_2d.num_faces())\n",
"print(\"Number of edges: \", mesh_2d.num_edges())\n",
"print(\"Number of vertices: \", mesh_2d.num_vertices())\n",
"\n",
"print(\"Coordinates: \\n\", mesh_2d.coordinates())"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 6,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "-cell 0 has faces : []\n",
+ " -face 0 has edges : [4 2 0]\n",
+ " -edge 4 has vertices : [1 5]\n",
+ " -edge 2 has vertices : [0 5]\n",
+ " -edge 0 has vertices : [0 1]\n",
+ "-cell 1 has faces : [1]\n",
+ " -face 1 has edges : [10 2 1]\n",
+ " -edge 10 has vertices : [4 5]\n",
+ " -edge 2 has vertices : [0 5]\n",
+ " -edge 1 has vertices : [0 4]\n",
+ "-cell 2 has faces : [2]\n",
+ " -face 2 has edges : [7 5 3]\n",
+ " -edge 7 has vertices : [2 6]\n",
+ " -edge 5 has vertices : [1 6]\n",
+ " -edge 3 has vertices : [1 2]\n",
+ "-cell 3 has faces : [3]\n",
+ " -face 3 has edges : [13 5 4]\n",
+ " -edge 13 has vertices : [5 6]\n",
+ " -edge 5 has vertices : [1 6]\n",
+ " -edge 4 has vertices : [1 5]\n",
+ "-cell 4 has faces : [4]\n",
+ " -face 4 has edges : [9 8 6]\n",
+ " -edge 9 has vertices : [3 7]\n",
+ " -edge 8 has vertices : [2 7]\n",
+ " -edge 6 has vertices : [2 3]\n",
+ "-cell 5 has faces : [5]\n",
+ " -face 5 has edges : [16 8 7]\n",
+ " -edge 16 has vertices : [6 7]\n",
+ " -edge 8 has vertices : [2 7]\n",
+ " -edge 7 has vertices : [2 6]\n",
+ "-cell 6 has faces : [6]\n",
+ " -face 6 has edges : [14 12 10]\n",
+ " -edge 14 has vertices : [5 9]\n",
+ " -edge 12 has vertices : [4 9]\n",
+ " -edge 10 has vertices : [4 5]\n",
+ "-cell 7 has faces : [7]\n",
+ " -face 7 has edges : [20 12 11]\n",
+ " -edge 20 has vertices : [8 9]\n",
+ " -edge 12 has vertices : [4 9]\n",
+ " -edge 11 has vertices : [4 8]\n",
+ "-cell 8 has faces : [8]\n",
+ " -face 8 has edges : [17 15 13]\n",
+ " -edge 17 has vertices : [ 6 10]\n",
+ " -edge 15 has vertices : [ 5 10]\n",
+ " -edge 13 has vertices : [5 6]\n",
+ "-cell 9 has faces : [9]\n",
+ " -face 9 has edges : [23 15 14]\n",
+ " -edge 23 has vertices : [ 9 10]\n",
+ " -edge 15 has vertices : [ 5 10]\n",
+ " -edge 14 has vertices : [5 9]\n",
+ "-cell 10 has faces : [10]\n",
+ " -face 10 has edges : [19 18 16]\n",
+ " -edge 19 has vertices : [ 7 11]\n",
+ " -edge 18 has vertices : [ 6 11]\n",
+ " -edge 16 has vertices : [6 7]\n",
+ "-cell 11 has faces : [11]\n",
+ " -face 11 has edges : [26 18 17]\n",
+ " -edge 26 has vertices : [10 11]\n",
+ " -edge 18 has vertices : [ 6 11]\n",
+ " -edge 17 has vertices : [ 6 10]\n",
+ "-cell 12 has faces : [12]\n",
+ " -face 12 has edges : [24 22 20]\n",
+ " -edge 24 has vertices : [ 9 13]\n",
+ " -edge 22 has vertices : [ 8 13]\n",
+ " -edge 20 has vertices : [8 9]\n",
+ "-cell 13 has faces : [13]\n",
+ " -face 13 has edges : [30 22 21]\n",
+ " -edge 30 has vertices : [12 13]\n",
+ " -edge 22 has vertices : [ 8 13]\n",
+ " -edge 21 has vertices : [ 8 12]\n",
+ "-cell 14 has faces : [14]\n",
+ " -face 14 has edges : [27 25 23]\n",
+ " -edge 27 has vertices : [10 14]\n",
+ " -edge 25 has vertices : [ 9 14]\n",
+ " -edge 23 has vertices : [ 9 10]\n",
+ "-cell 15 has faces : [15]\n",
+ " -face 15 has edges : [31 25 24]\n",
+ " -edge 31 has vertices : [13 14]\n",
+ " -edge 25 has vertices : [ 9 14]\n",
+ " -edge 24 has vertices : [ 9 13]\n",
+ "-cell 16 has faces : [16]\n",
+ " -face 16 has edges : [29 28 26]\n",
+ " -edge 29 has vertices : [11 15]\n",
+ " -edge 28 has vertices : [10 15]\n",
+ " -edge 26 has vertices : [10 11]\n",
+ "-cell 17 has faces : [17]\n",
+ " -face 17 has edges : [32 28 27]\n",
+ " -edge 32 has vertices : [14 15]\n",
+ " -edge 28 has vertices : [10 15]\n",
+ " -edge 27 has vertices : [10 14]\n"
+ ]
+ }
+ ],
"source": [
"for cell in cells(mesh_2d):\n",
" print(\"-cell\", cell.index(), \"has faces :\", cell.entities(2)) # cell and faces are identical in 2d\n",
" # try this on a 3d mesh \n",
" for face in faces(cell):\n",
" print(\" -face\", face.index(), \"has edges :\", face.entities(1)) \n",
" for edge in edges(face):\n",
" print(\" -edge\", edge.index(), \"has vertices :\", edge.entities(0)) "
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 7,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ " \n",
+ " FEniCS/DOLFIN X3DOM plot \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"mesh_3d = BoxMesh(Point(0,0,0), Point(1,1,1), 5, 5, 5)\n",
"mesh_3d"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# Save mesh in VTK format\n",
"file = File(\"mesh_3d.pvd\");\n",
"file << mesh_3d;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Besides those simple geometries, we will use meshes created from medical images:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 9,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[,\n",
+ " ]"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAAD4CAYAAAAw2nCGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXrklEQVR4nO3de5BV5bnn8e/TN6DB5hIu4WZgJjQVxcuMrZg6JDAmnkBCnRwpSaFxojkqw81SxyACSh8knQAK48SJVpwca7TiiWMSSYwDOeIFC8sJtxQeQUVRc8AIiOnIACYNTT/zR+8mTbv7utfa79pr/z5VFLvX2nutZ9v+WGuv/a7nNXdHRNKnJHQBIhIPhVskpRRukZRSuEVSSuEWSamy0AV01eDBg33MmDGhyxBJlB07dnzo7kOyrSuYcI8ZM4bt27eHLkMkUczs39pbp9NykZRSuEVSSuEWSSmFWySlFG6RlFK4RVJK4RZJqYL5nlui8etf/5rf/e53ZyybPHkyU6ZM6fE2ly9f/oll06dP5+mnn+7wdbW1tT3ep3TOCuV+7pqaGtcglu45ePAg27dvZ9q0aWzYsIGGhgZ27dqV9bnnnnsuV155ZZe2++qrr/Lkk09GWeppQ4YMYd68ebFsO43MbIe712RbpyN3gTtw4AAff/wxx48fB8DdMTMA1q1bB8Abb7xxen17du/ezaWXXsrIkSNPvz6b/fv3xxZsgMOHD7NlyxZ69+5NSUkJ5eXlNDY20tjYSGlpKRUVFYwePZrKysrYakgLHbkLRMv/3E1NTbT8zo4fP859990X6X5uueUW+vfvn3XdsWPHWLNmTaT766lly5Z1+I9QsdCRu8CdOHGC73//+7Hvp7y8nKqqqqzrjhw5Evk/JLm4++67FfBOKNwJlO0CVT4sXryYhoYGVq1axZQpU5g8eTLQfMr+85//PEhNHbn77rs/sUwX6f5Kp+UJcfDgQX70ox+FLuMMl112GWVlZTzzzDOhS+m2JUuWUF5eHrqM2HV0Wq5wB7Rp0yYOHjzI4cOHqa+vD11O6owfP54vfvGLjBgxInQpsdFn7oRwd/785z8DzVe5X3zxxcAVpduePXvYs2cPS5Ys4eTJk0V3hV3hjlnrM6NsnxElft/73vcAWLRoEb169QIoigtxCneMtm7dyoYNG0KXIRmrVq064+e0X3yLPdxmVgpsB/7g7tPNbBDwv4ExwO+Bb7j7n+KuI59Wr159+vRbkmv58uWpDng+bhy5GXi91c93AM+5+zjguczPqbF8+XIFu4CsXLmSpqam0GXEItYjt5mNAr4G1AH/NbP468CUzONHgE3AojjriFtDQwNvvfUWO3fuDF2KdFNDQwMbNmxg4MCBnHXWWQwbNoyhQ4eGLisScZ+W3wfcDpzVatkwdz8A4O4HzKxg/0u2DAVduXJl6FIkB22/Yl26dCllZYV/OSq203Izmw584O47ctjGbDPbbmbbDx8+HGF1uWtqamLFihV897vfDV2KRKyuro4TJ06ELiNnsQ1iMbPvA/8ZaAR6A1XAk8DFwJTMUXs4sMndx3e2vSQNYmlsbKSuri50GZIHF110EdOnTw9dRrs6GsQS25Hb3Re7+yh3HwPMAp5392uAp4BrM0+7FvhVXDXE4dChQwp2EdmxYwe//OUvQ5fRIyE+WKwEnjCz64F9wMwANXTba6+9xpYtW9i3b1/oUiTPXnnlFQYPHsykSZNCl9IteQm3u2+i+ao47v5H4Ev52G9U6uvr+dnPfha6DAnoueee44ILLqCkpIS+ffuGLqdLCv+SYMyOHj3K/fffH7oMSYC1a9cCUF1dzVVXXRW4ms6p+2kH9u3bd/oXKtLizTffLIhvSXTkzmLjxo28/PLLocuQBDt16lTih6/qyJ2Fgi1dleSPbAp3G2+++WboEqSA1NfXs23bNk6dOhW6lE9QuFs5cuQIP/3pT0OXIQVm/fr17NjR44GYsVG4M/70pz8lqrunFJYNGzYk7hRd4aZ5Bo0f/OAHocuQAldfXx+sc202RR/ubdu2xTqDhhSfpHxNVvThXr9+fegSJGVOnTpFEroKF3W409qBQ8I7duxY6BKKO9wrVqwIXYKk1Nq1azl69GjQGoo63CJxWrt2LZs3bw62/6INd5Kuakp67d27N9i+iy7cjY2NCrbkzb59+9i1a1eQfRdduA8dOhS6BCkyv/jFL4LsN9Zwm9loM3vBzF43s91mdnNm+SAz22hmb2X+HhhnHa2dPHkyX7sSOS3EV2NxH7kbgdvc/XPApcB8MzuHgBMTHDlyJF+7EjktxDxxsd7PnelP3tKj/KiZvQ6MJNDEBJs3b+b555+PezciWbVc65k7d25eJj7I22duMxsD/AdgC20mJgCyvtOo+5a/8MILOW9DJFf5uu6Tl3CbWT/gF8At7v7/uvo6d3/I3WvcvWbIkCE51fD2228nYkigSL4Gt8QebjMrpznYj7l7yx0ahzITEpD5+4O46/jJT34S9y5EuuSjjz7Ky9DnuK+WG/BPwOvu3rrTYF4nJtARW5Jk27ZtrF69Ovb9xN0g8W9onlLoVTNrmQJzCXmcmMDdg1ypFOlIQ0ND7PuI+2r5S4C1szovExM0nzyIJM8Pf/hD5s+fH9v2Uz9CTUNNJak+/PBDHn744di2n+pw5+PURyQX+/fvj+2aUKrDfeDAgdAliHQqroNQasO9detWHnnkkdBliHRq1apVfPzxx5FvN7Xh1lFbCkl5eXnk20xluLds2cLOnTs7f6JIQsQxqCWV4f7Nb34TugSRboljOqLUhfsvf/lL6BJEEiF14V61alXoEkS6TaflIim1Zs0aXnrppUi3qXCLJMRzzz0Xab81hVskQaLslKpwi6SUwi2SMFGNNVe4RRImqv4DwcJtZlPNbI+Z7TWzvLU2FikWQcJtZqXAD4FpwDnAVZl+5iISkVBH7kuAve7+jrufAB6nuZd5TkJN2yKSRKHCPRLY3+rn9zLLztCdvuVNTU3BJlwTSaJQ4c7W2OwTlwij7FsuUmxChfs9YHSrn0cB7weqRSSVQoV7GzDOzMaaWQUwi+Ze5j2mLqciZwoSbndvBBYA/wK8Djzh7rtz3GYUpYkkQhTfdcc9KUG73H09sD6ibbFu3booNiWSCO5OY2MjZWU9j2gqRqj98Y9/1JVykTZSEe58TKomUmhSEe5evXqFLkEkciUlucUzFeEuLS0NXYJI4qQi3JWVlfoqTFIn14+bqQh3fX29vgoTaSMV4RaRT0pFuHP5LlAkrVIR7oqKitAliEQu1wvFqQj3yZMnQ5cgErlcpxhKRbjjmGdJpNClItz6Gkzkk1IR7gEDBnDeeeeFLkMkUVIRbjPjiiuuCF2GSKTq6uo4dOhQj1+finCLpNXRo0d7/NrYwm1m95jZG2b2r2a2zswGtFq3ONOvfI+ZfSWK/WmEmqRRIsMNbAQmuPv5wJvAYoBMf/JZwLnAVOCBTB9zEWkjl4vFsYXb3Z/JtFMC+C3NTRChuT/54+7e4O7vAntp7mMuIm3kckaar8/c/wBsyDzuUs9yEcntnu6cBmWb2bPAp7OsWuruv8o8ZynQCDzW8rIsz8/6z5OZzQZmA5x99tkd1pLrje0iSVReXt7j1+aUCHf/srtPyPKnJdjXAtOBb/pfzy+63LO8u5MS1NbW5vJ2RBKnsbGx8ye1I86r5VOBRcDfufvHrVY9Bcwys15mNhYYB2yNqw6RQpbLZ+4475X8H0AvYGPmit9v3X2Ou+82syeA12g+XZ/v7hocLpJFnz59evza2MLt7p/tYF0dUBfXvkXSIpebonQVSiTBcjlyK9wiCVZZWdnj16o/kUhC3XbbbfTr16/Hr9eRWyShevfundPrFW6RlFK4RRLIzHLu6qtwiyTQsmXLct6Gwi2SUgq3SEop3CIJdOzYsZy3oXCLJNCaNWtyuiMMUhjuKVOmhC5BJBK59uNPXbgnT57MZz/b7j0rIgUj1wYkqQs3wDe/+U2mTp0augyRnGiusHZMnDiRW265JXQZIsGkNtwA/fv35/LLLw9dhkgQsYfbzL5jZm5mg1sti3xSgvZ8/vOfj3PzIrHo1atXzsNPY73l08xGA5cD+1otaz0pwQjgWTOrjqvVkmYAlUJ0xx135LyNuI/c/w24nTNbF2tSApE8iLP76d8Bf3D3V9qs6vKkBGY228y2m9n2w4cP96iOF198sUevEyl0sU1KACwB/jbby7Isy9q/1d0fAh4CqKmp6XaP14aGBjZt2tTdl4mkQk7hdvcvZ1tuZucBY4FXMp95RwG/M7NL6MakBCLSc7Gclrv7q+4+1N3HuPsYmgP9H939IHmclEBTDEkhGjZsWCTbyXuDxHxOSlBaqpmBpfDMmTMnku3k5dCWOYJ/2OrnOnf/9+4+3t03dPTaXJSUlETS0UKkEKX+vNXMuOuuu0KXIZJ3qQ836LO3FCf9Xy+SUgq3SIJceOGFkW1L4RZJkJ07d7Jx48ZItlU0c4XV1tbS1NTEgQMH+PGPfxy6HJF2vfzyy/Tr1y/nOxqL6shdUlLCiBEjmD9/fuhSRDr0zDPPqEFid5kZAwcODF2GSIe++tWvajqhntA93pJ0F198cc7bULhFEmbQoEGRbKdoLqi11tTUFLoEkaxqa2sj21ZRHrlFkmjAgAGRbq8ow11aWsqsWbNClyFyhqiHSRdluAEef/zx0CWInCHXSQjaKspw6zO3JNGnP52tY1nPxRpuM7sp05t8t5mtbrU8b33L22psbGTFihX53KVIl4wdOzbS7cV2tdzM/hPNbYzPd/cGMxuaWZ7XvuVtffTRR/nYjUi3XH/99YwaNSrSbcZ55J4LrHT3BgB3/yCzPGjf8sGDB3Prrbfma3ciXRLVd9utxRnuauALZrbFzF40s5YhN3ntW55NVVUVM2bM4PLLL49kJJBIEsXZt7wMGAhcClwMPGFm/4489i3vyHnnnXf6cX19PW+//XaUmxfplsrKysi3GUvfcgAzmws86e4ObDWzJmAwCetbvnnzZgVbgopyVFprcZ6W/xK4DMDMqoEK4EPy2Le8K55//vlQuxYBoPn4F704x5Y/DDxsZruAE8C1maN43vqWd+aDDz7o/EkiBSq2cLv7CeCadtbVAXVx7burHnzwwdAliMSmKEeotdD4ckmzog73+PHj+cY3vhG6DClycfUXKOpwA3zuc5/j29/+dugypEjFdaUcFG4Azj77bG6//fbQZYhESuHO6NOnD9/5zndClyESGYW7lb59+/K1r30tdBlSJL71rW/Fun2Fu42amprQJUiRiPoWz7YUbpGUUrhFAhg5MuuNkJFSuLOI8+sJkerqam644YbY96Nwt6N3796hS5CUmjBhQl72o3C347rrrgtdgqRUdXV1XvajcLdj2LBh3HjjjaHLkJS588476dWrV172pXB3oF+/fqFLkJTJ5zx1CncHqqqqWLp0aegyJEWinnigI7GF28wuNLPfmtnOTJPDS1qtC9a3vLvKyspYtGhR6DIkBRYuXEh5eXne9hfnkXs1sNzdLwSWZX5u27d8KvCAmZXGWEfOevfuzdVXX80555wTuhQpYHE0QexInOF2oCrzuD9/bYIYtG95T40bN46ZM2dy7rnnhi5FCtC4cePyvs84w30LcI+Z7QfuBRZnlne5b3kSXXDBBaFLkAJ09dVX532fOYXbzJ41s11Z/nyd5hlHbnX30cCtwD+1vCzLprK2f4xrUoJcjBs3Tp/BpSDE2bf8UeDmzI8/A36cedzlvuVxTkqQC80SKt0RajhznKfl7wOTM48vA97KPE5U33KROPXt2zfYvuPsW34j8N/NrAz4CzAbwN0T07dcJG5xTPDXVXH2LX8JuKiddYnoW95Tffr0CV2CFIjLLrss2L41Qq0HzIxly5aFLkMSbuHChYwZMybY/hXuHlLApTNxzQHWVQp3DhRw6YjCXeDyeZePFI6LL744+F2FCncE4m5RK4Vn2rRpoUtQuKMwduxYrrrqqtBlSIIk4YxO4Y6Ivh6TFklpsKlwR2T06NEsXry48yeK5InCHaGKigrmzp0bugwJKEk3FSncERs6dGjoEiSQ2bNnJ6oltsIdg0996lOhS5A8mzZtGsOHDw9dxhkU7hgsWLCAYcOGhS5D8mjEiBGhS/gEhTsmc+bMYeLEiaHLkDy49dZbGTVqVOgyPkHhjtHUqVO56aabQpchMbr55pupqqrq/IkBKNwxy9fsEpJ/V155JQMGDAhdRrsU7piF7MQh8Up6J9xcGyTONLPdZtZkZjVt1mWdeMDMLjKzVzPrfmBJGKcXs6SMWJLoTJo0KXQJncq1E8suYAbwo9YL20w8MAJ41syqM+2UHqS55dJvgfU0T0ywIcc6Eq+2tpbly5eHLkNydNNNNwVtndQdOR253f11d9+TZVXWiQfMbDhQ5e7/15tvdn0U+PtcaigktbW1lJYmenIV6cBtt91WMMGG+D5ztzfxwMjM47bLs0pi3/JcLVmyhC984Quhy5AeyPd0QLnqNNydTDzQ7suyLPMOlmfl7g+5e4271wwZMqSzUgtCSUkJkyZN0iCXAjNx4kRKSgrr+nOn1br7l919QpY/v+rgZe1NPPBe5nHb5UWloqKCOXPmhC5Duqi6upqpU6eGLqPb4mpt/BTwz2a2luYLauOAre5+ysyOmtmlwBbgW8D9MdWQeLrIlmzXXXcdn/nMZ0KX0WM5hdvMrqA5nEOA/2NmO939K51MPDAX+F9AH5qvkqf+SnlHWr4me/rpp9mxY0fgaiRNX1ta6A6NXVVTU+Pbt28PXUas6uvruf/+oj2RCWbGjBmMHz+eioqK0KV0m5ntcPeabOvinE5IumnQoEFMmzaN48eP09DQwJYtW0KXlGrnn38+I0aMYMKECYnoeRY1hTthLrnkktOPx48fz6OPPhqwmvSaNGkSX/rSl0KXESuFO8HGjh3LvHnzeOCBB0KXkhqzZ89OXFOFuCjcCTdkyBCWLl1KXV3BzpuYCMuWLUvlqXdHFO4CUFZWxl133YWZ4e489thjvPPOO6HLKhjFGGxQuAtGy+goM+Oaa67h6NGjuDvl5eXcc889gatrNnv2bPr27Yu7Y2anp9NZsWJFXva/cOFC3J3Gxkag+b/VWWedVZTBBt3PXZDMjKqqKvr3709lZSUzZ84k5PDcQYMGMX36dIYPH366rqqqKkpKSigpKWHBggWx13DjjTdSWVlJ37596d+//+kaijXYoO+5U+/dd9+N9Ir7hAkT6N27Ny2/i+uvv75L/cNOnDjBfffdx4wZM3jhhRd4//3oRh3Pmzcv6D9uIXX0PbfCLQDce++9HD9+vN31N9xwAyNHtnsDX85aD8O98847dWtsFyncIinVUbj1mVskpRRukZRSuEVSSuEWSSmFWySlFG6RlFK4RVJK4RZJqYIZxGJmh4F/C7T7wcCHgfYdQrG9Xyjc9/wZd8869rZgwh2SmW1vbxRQGhXb+4V0vmedlouklMItklIKd9c8FLqAPCu29wspfM/6zC2SUjpyi6SUwi2SUgp3O8zsH83sD2a2M/Pnq63WLTazvWa2x8y+ErLOqJnZ1Mz72mtmd4SuJy5m9nszezXzu92eWTbIzDaa2VuZvweGrjMX+szdDjP7R+CYu9/bZvk5wE+BS2iewfRZoLrVRIcFy8xKgTeBy2mebnkbcJW7vxa0sBiY2e+BGnf/sNWy1UC9u6/M/MM20N0XhaoxVzpyd9/XgcfdvcHd3wX20hz0NLgE2Ovu77j7CeBxmt9vsfg68Ejm8SPA3wesJWcKd8cWmNm/mtnDrU7RRgL7Wz3nvcyyNEjze2vLgWfMbIeZzc4sG+buBwAyfw8NVl0EinpSAjN7Fvh0llVLgQeBFTT/T7ACWAP8A5CtEXZaPtuk+b219Tfu/r6ZDQU2mtkboQuKWlGH292/3JXnmdn/BJ7O/PgeMLrV6lFAdE24w0rzezuDu7+f+fsDM1tH80eSQ2Y23N0PmNlw4IOgReZIp+XtyPxyW1wB7Mo8fgqYZWa9zGwsMA7Ymu/6YrINGGdmY82sAphF8/tNFTPra2ZntTwG/pbm3+9TwLWZp10L/CpMhdEo6iN3J1ab2YU0n5b+HvgvAO6+28yeAF4DGoH5abhSDuDujWa2APgXoBR42N13By4rDsOAdZmphsqAf3b335jZNuAJM7se2AfMDFhjzvRVmEhK6bRcJKUUbpGUUrhFUkrhFkkphVskpRRukZRSuEVS6v8DVpqWOV8cBbMAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
- "path_to_file = \"../data/domain_brain_2d.h5\"\n",
+ "path_to_file = \"data/domain_brain_2d.h5\"\n",
"# load mesh from HDF5\n",
"mesh_brain_2d = Mesh()\n",
"hdf = HDF5File(mesh_brain_2d.mpi_comm(), path_to_file, \"r\")\n",
"hdf.read(mesh_brain_2d, \"/mesh\", False)\n",
"\n",
"# plotting via ipython 3d rendering (slow)\n",
"#mesh_brain_2d\n",
"# or via matplotlib in 2d\n",
"plot(mesh_brain_2d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A mesh can be devided into subdomains:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 10,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAAD4CAYAAAAw2nCGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbtUlEQVR4nO2dfewmV1XHP2dbQaMUynYLpS1uja0JImnCj7opMVRgaUV0wXRNMWiVxtaGyosxQOEfDCE0iPiGkq2hEROwFAS3kWLZotWQdWl3SbUttXSFvjd0WxA1Rpptj388z2xnZ+fl3pl7Z+7cOZ9ks79nnnlm7szc7z33njn3XFFVDMPIj01TF8AwjDiYuA0jU0zchpEpJm7DyBQTt2FkyvFTF8CVk046Sbdu3Tp1MQwjKQ4cOPCYqm6p+2424t66dSv79++fuhiGkRQicl/Td9YtN4xMMXEbRqaYuA0jU0zchpEpJm7DyBQTt2FkionbMDLFxG0YmTKbIBYjHNs37Tzq856nPhP0eMUx67aHOqfRjYk7Y7Zv2nlERFWhPX7ZubX7uR536D5N35vow2Hd8gVSFraRLzKXNEsbGxtqseXH4mJFu3j8snPZvGuvk9UMcb6hmHV/GhE5oKobdd9Zt3wmxBTVnIQN/sOIpWLiTpC+Iip3tzfv2uv8u7JQinPXbRuLtmFDcV1NTjzjaUzcCeEiItfxcnW/NrEXlnD7pp1Hflfe1peiuz8WdQ3TkjGHWiKk0uU18sEs94T4Cjqml7soS9nSNpXPpxwuPYgQFr58jHK5l2zFzVs+IkOtcwxxN4mqSXAxG5jq+R6/7FyO33GIw7u31O5TV5am68lV5JN6y0XkOGA/8JCqvk5Engt8GtgK3Av8sqp+N3Y5piB0V3vzrr29xdX0Wx+r2efcx+84BHCUQNuOXy7L5l17eRy/czY5FUNH5c2BMcbcbwPuKn1+N/BlVT0T+PL6c1Zs37TT2TnmKxjf7uvmXXs7f9NVhj7lrHL8jkPH/PMpS/U6XMrTts8SfBxRLbeInAb8PPAB4HfWm3cA563//gRwM/CumOUYk7ZK01TZ+r7CcqGvd93ltz5Wue33BU3HGXJP2u5t7u/LY1vuPwLeCTxV2vY8VX0EYP3/yZHLYBiLJJrlFpHXAY+q6gEROa/nMS4FLgV44QtfGLB0YQnhVZ4bZatbdXoNOebh3Vui9WTqPPc5j8WjectF5IPArwKHgR8ETgA+B7wMOE9VHxGRU4CbVfUnuo6Xorc8lqhjO7hCCbFK3XGbxtZNVI8xpse+ONecBN7mLY/WLVfVK1X1NFXdClwE/IOqvgm4Hrh4vdvFwO5YZYjJHK11ITRfwaVEzPtbjs5zdYqmzBQRalcB20XkHmD7+vOsqHvobR7lLg9x3bFi41OeOkJY/zGOWcblOVTv/ZwFPkqEmqrezMorjqo+DrxqjPPGwPdhz9lKzoVQQ43csPBTR/p0w11f9bgcd6hjaYxGpu0c5Wtv2q+PSKtDjWo02/duONP7XNVgmrm+MjNxOxBC2FPgKrbYxDpX3fUdv+MQz37tPUedu09jUvXYz3HGmYm7haaxdR9CVvBqGaZqSIpr6jr/2OWr3p++Ai8frzopZQ4iN3E3EFLYTfiEU7Y563wJbUnbxJMKQ8tY3P+yyFMXuIm7htAeUpd3t2ORsuMpdtlcxv1dlEWeusAtWUNkUhbT3Eg5hiBFzHJX8O2ON02eiDEJouncrozp2PLFZdw79BhNDHE8pmy9zXKXqEvc7+oRdwkKaUuM0HbsUIQ+5tBAmC6K+xVD2K77N11jdVGHFINdzHLTz3lWfuBDLOJYwk6ZNp9E13OYctiT+vtwE3cNPp5p18rVZwJEUwPSV/ypjv+rr6NSjx8ok7LAF98tn6o7Zc4hf1ITduqY5a7gK7o+IZN9LFOIij2k5xGbtvs4tahdnX1TvuKsw8Q9Ai4PfeoK3AfXCDVX+hxnaKqnLnznA6T0/nvR3XLfLnndg257+HXJCXPsjh/eveXIv9zwbXDKc8KnZrHibrr5XVa2WonnaHHLhC5/jiJPcT6+C4vslodoVcuxyj7j7qbXPKG7uHVMMRMspcYvht+i7Z5OPckkquUWkdNF5B9F5C4RuVNE3rbe/lwR2SMi96z/PzFmOQxjiURdTmidAPEUVf2aiDwLOAC8Hvh14DuqepWIvBs4UVVbc5eHTJBYF4nWl6G5t8eavjllV3lM6z22x73L5wJxLXdbgsRR1woTkd3AR9f/vDKghhJ3SGFDc6XxEXlMgacw/k2paz6UPvdz86690QSehLhFZCvwz8CLgftV9Tml776rqsd0zSt5y1963333DS5HWdwhHB++4i7TFmY5xmSMKZij0PsKuo7QIp8ktXGlAD8C/A3wdlX9L9ffqerVqrqhqhtbtgyvqGO+npi6Ert0F41ufIRdvPpM5f5GF7eI/AArYX9SVT+33vztdXe8GJc/GrscxvRM3eD54ivs1Ii9EKAAHwfuUtWPlL4qFia4ihkvTNCFT/RUdXpjrNRFU1TCuYkaup9ZimKuEvs998tZLSl0u4jctt72Hlaivk5ELgHuB6L3l4c40nwqZwrj3LYyTBFgMYcca2WGDmnaliEe8913VHGr6lcAafh6koUJXCt3zCwnc2j1Q5NiUIsvPsk2mvYb89kvIvy0aC19rNbYFjjEAvdzIIWeTR+GCrvKGM7d7MNPh9zEGNbGNWnDnC1cF3PrpkO44cyY2VMXYblTZAlW2piWrC13qK6Pi6Vx8YxXPeIhFxowwhEqP97ULMpyT+nIcs3kGWrmUuzMpEsh1j0cY953tpZ7qsnyTdM/xxB113FTsUJza3SanmkqCy02ka24m/DJhT0WY1X2oULfvGtvUvdtTHyfke8CgzHIUtxTprjpu770mPiWsVwB+wp8btY6BLHzu3WxqDF3wVKtj7EssrTcVfqkK+6LT8qlpVizpVxnE33SXwc57+hnHJkxhV0+hq8DJmfmGLTSlz4ijhXMkr24fQhZAXOtzDakOZpQFjnGhJKsxT2nNadSxgRdj8/ssSnuYbYONRO2sXSyttxNLGkMaMShbLVd3lXXWfLYEZOTWW4RuUBE7haRg+v0xqOS48oYrozdsOXUkFbrTV+BNv0uZIzGJOIWkeOAPwN+DngR8EYReVHIc7je9CULfAzR5SbsKin7I6bqlp8DHFTVbwKIyLXADuDrQw9cbvlcQ01zyBLSl9DvYJd4D6v1K5VMO1OJ+1TggdLnB4Gfru5UyVveedCuxf1SbmWnJJTAcxe2633qU89iNAhTjbnr8qodszpC6LzlbeReMY3lMZXlfhA4vfT5NODhicpiwma6EMm54bPKpyuxuvFTiftW4EwROQN4CLgI+JWYJ7QueXyW+Iox5Tndk4hbVQ+LyBXAjcBxwDWqemfIc8RKYZwzZr3D4mrlY73zniyIRVVvAG6IcWwTdj9M2HFpS5YRIytqNuGnUyZoMIwUyUbcBWa1DWNFduJuG7tYVtB27L6MR9O9Dun4zWbiyJ6nPnOka16X58sqbnzsHq9IxXeRjbjLxHrt1fTQcqrUUyf1mzsp3bdsxN212F/xDrbr5jcJNaWHNga+r8VyauD6klodyW7MbRjGimwstwsuLesc8o4b86DNuTvGTLJFiTs0JmwD6g1Cl1jHmBZq3XIjKVIbt86Z7MS9edfe2laxaXtfzGofTUhRzikFVh+rPRbZidswjBXZjLnLQSyQTus5V/pYzupvfHo3+87+LNtuu/CY41kPqT/ZiNsIwxTd4X1nf/bI/9u4cDZdcqiPB6iLtehapMC85T0JGq+bsSVxFVWodboLUVe3beNCr/LMga77VZ7yCWHWD4smbhH5feAXgCeA/wB+Q1X/c/3dlcAlwJPAW1X1xtDnjxGCmrOwXbHhzrG0zdP2JWTihpgOtT3Ai1X1JcA3gCsB1vnJLwJ+ErgA+PN1HvOkMWEfLWwTeT11sw59ZyKGMkzRxK2qX1LVw+uP+1glQYRVfvJrVfX7qvot4CCrPOaGcYRcuuRTGoWxxtxvBj69/vtUVmIveHC9LTnmbK3rvM9j0nXv6sbb1d/PWeB111+3LeY1DhK3iNwEPL/mq/eq6u71Pu8FDgOfLH5Ws/8xOcvXv/ValKBwQGzftPOobqNvN2eOoj6mkpztL3AfQTXd0zneuymJ2YgNEreqvrrtexG5GHgd8CpVLQTsnLNcVa8GrgbY2NiobQCMdgoL6SrytvncXY1kDGGn/q475Tn+Mb3lFwDvAl6hqv9b+up64FMi8hHgBcCZwC2xyuFLn4dSfsApPNQ6+lhxcO82+gaszJ2u+5JCnYg55v4o8Exgj4gA7FPV31LVO0XkOlaL/h0G3qKqT0YshzFDUh5vp1y2MtHErao/3vLdB4APxDp3GZ/xdlcL6zofPFXr3YfJrE5pLJrSPQ0RllsQ+5oWEaHWRQhRz4GpPehzx6UeNL3/rzMyRaMVy6mWrbiHOn+6bvYQb/yUTC3wOY+3u0TYFtjTFLKb7KuwlOmb3rjpZltEltFFV+ho3XdNk0xCLCeUrbj7UBa2i5jnZLHLlK2nddP96Iojr04A6SKm0cg2WYOv8HIZV+eCPY/hLN5yVyvR0JY0Fa/u3JlD+KlrwM9UQ7pFi7uPsC3s0qjSFfDjO40zxHgbMha3T0KBrhs/RdjlWIzlPZ+zl9wVl/W3m7BMLAPxyVQ596WAt9124WwFNSQXWyqETODQF3OoGUamLMpyV3F971hH6tbk8O4tcLbbvrG75j49iBxfzblYcd9XaE7nDXakhCjSHLcJ1TfCLHUxG/PAReghkiNCpuKGY/OY92kZTdBGlZD+gNiv+7IVdxVfx1kOwp6zU63M1LPC2gQY2vkXMvtp1uKuWu8quTvdfMbdxtH0tag+SRpie9Gz9ZYX7HnqM8GCAgxjTkQXt4j8roioiJxU2naliBwUkbtF5PzYZaijzWqH7gKO8Z4zh2FEKvisvNLWhS5WK53qPXfUbrmInA5sB+4vbSsvSvAC4CYROSt2qqWuLnrB0BxqbfssVYAhxv7FPY59D33m8Tdta5u37VL+uYSf/iHwTmB3aduRRQmAb4lIsSjBv8QsSCHsUBZ7SLqdFEU+dRKHFOgj7Lb9+iRn2LxrL3zsHU7n6SJm9tNfBB5S1X9dJ0gscF6UwDdveRNdwnYVW+qzlOooBOtiOW2ed1jaRD4G0RYlAN4DvKbuZzXbanOSh8hb7mKxDSNHoixKICI/BZwBFFb7NOBrInIOHosSpESf3FkpNChHyj3xKzGfHsRUdAWV9J2jPVW+vSjdclW9HTi5+Cwi9wIbqvqYiCS9KEGZvsnw5ox1x91oEmlXvRiz3owexDL2ogTl9cPA3ZPpMz20jhSsdoqkHjXnu9JKlaERZiFjMkYRt6purXwebVGCKi5JHIYI20SdB0NEXq0DlmYpMi7vuYekNTZRu5O69S7T9ialz4qoYwo9+/DTKr4izGFsneJ7dSM+i7HcfcilK+46gWQOzrSUIv3qytFlzce04osTt0uoYC6iLpNSVziVcsTAZ6wecnpnHYvqltd5Iqs3d87CLhaVq6Orso1ptefQQxhK8Szangms6tHjl50bZfbi4ix3nWMtxiuuOSTVn5KUehJTs3nXXraX6uBcJo4kiesMsSpDJ51MPVZsEtRUlrRP1NrU9zAG1XoVKofaorrlhrEkFmm5fem74khqVubw7i1s48LkZn+Vy5DjUMb3mqxbHpk5pztuG+8f3r2Fjd2XH7Xv1OQoaPDL6FIwu/DTOZHzumAw/RzjOnJ0PvoKO0aePxO3I3MTdZNgClGndj1dAk+tvKGI2cgu1qFmWVGN3FmsuLdv2nnM67CcFygoE+t6hnatXQM/5oBP+WNFqVm3vIUlVTBfqkJeVdB7eJynG8ih58/p/g+Zcdj7/NGODIjIbwNXsErK8AVVfed6+5XAJcCTwFtV9caY5ajimlfNJVtpShMZxqBrrnvZYeez+kbuTOE0jJn99GdZpTF+iap+X0ROXm+fJG95Ha5L+DaJPDcPbxN9rE41KUbKaZ2nJpbvJ6blvhy4ap2fHFV9dL19krzlZco3s08YKiynsk7RncyRoWm7+hDToXYW8DMi8lUR+ScRedl6+6nAA6X9WvOWi8h+Edl/6FAcEZnX3MiVmHnLjwdOBLYBLwOuE5EfY+S85bHI3WLHYmk+Cphu+BYlbzmAiFwOfE5VFbhFRJ4CTmKmecsNwwfXpYli9hpjjrn/FnglcLOInAU8A3gMSC5ved9xdxu5RFv19fK2vYnIzXqn6liNKe5rgGtE5A7gCeDitRUfNW95F6GXG3Jd8bPMHCt67BRBc8FX2FksSqCqTwBvavhusrzldYwpbMMYi8VHqFWzsvTtTg4Rdm7dVHBb/CFlXJ5n6s9ssbHlZQqnRt/K6FIRNu/aO9tubA6x3j64NtR9fRFjNXombk8O795y1EN1FXaOjFlRx2LMoVXs+IrFd8sL6pImdnXRm+gKy3Q5/twZmkwyV2w5oYmoRqv1eRB91202jNCYuGsYOgY3jBQwcbcwZydYquTUJfetG7ESITZhY+5A5NIIhHIozf1VWBd1c9ghrd6eiXtCUqoIMH0QTgrTaF3DbYtn17TWXPXZTtH4m7gbKHvPu6xQjLXGcqUuYGfqRqVKiKwpTXVizOnFNuY2JiU1YeeEWe4W+i4YWEfq733HFFnXuVLpnpfL0kTKE2jMcndQdKOaPOcuD3aJwg4x/JjSqlcjEdtwjdQbO+OPWW5Pqg6TtpY79zzoqVqsofRtVFKz4ibunrQ9xNQtdZUmB1L1Gps8xEaaROuWi8jZIrJPRG5bJzk8p/TdlSJyUETuFpHzY5UhFKGSKKYo7CaahiCuwrYG4GimSMIZ03J/CPg9Vf2iiLx2/fm8lPKW+1I8oBhpmQwjNDHFrcAJ67+fzdNJECfPWx4Tl5lkLp7Yqa18auPHVJjTPYkp7rcDN4rIh1l1/4tafyqwr7RfY97yVPGdHlolxXe7dePuUIE7c8PVB+HKVHnxB425ReQmEbmj5t8OViuOvENVTwfeAXy8+FnNoWpzko+xKEEoQk8ymaIB8Fllc4xIu6nfc1fPP7foQlklJI1wYJHvAc9RVRURAb6nqiesFwFEVT+43u9G4H2q2tot39jY0P3790cpax+6xt0hKsLYlbtvg+Lz5sB3QsnUw5O+STlgHIstIgdUdaPuu5hBLA8Dr1j//UrgnvXf1wMXicgzReQMEshbboxDzl35FIk55v5N4I9F5Hjg/4BLAVQ1qbzlxoohw4A259tQQU+dGbZtEkkKM7/aiJm3/CvASxu+SypveR9yei025vi+z3BliMCHrhGeovPTFYtQG0jIySVlxrJWY1TesR1RTdfkMyFl6H1JYeVYE7dxFE0hp3PBZzmnUAtMpNYdLzBxByC09a6rdHUVLrR1d5n1Vif2vnH2rsTqxYSYWpqqsMGmfM6COY/7jOkwyx2IqvUOlSCwTdhTJDVI2VKFaARDHCOF8TaY5U6OslB91qzqWynL50txfD1Gw1W+f1MHzYTExB2RIiS1T2hqk1hjW84iq0gooQ8pr6/QfPZvK9cQgaditSFi+GloUgs/baLNseaTmcUneULbcXwYe1HDrgYk1nvptjzjPplZ6+7F2OKeKvx0kYRK7LB0Uugez92RaeJOgGolCjnd0AcXQY2xbO8QYYfqvfQRdmqNunnLI+H77jsVh45PSt8yKXvR26h7q+HbtU9N1AVmuSNSfehDBNAVPBK6UXCd192XmNbft9y5Lvho4h6ZLu9533FezPFhDJFPKaamhqWrTHXPLlWrDSbu6IR8+Cm+h06Rrvf+fRqWFDzjvpi4DSfm7jleIibuEWhq4Yd2Tcey5KGF3eZxH+NdfRt1Q6Y5Wm0YniBxp4jcKSJPichG5bvahQdE5KUicvv6uz9Z51fLnqb33z4C79p3SBhq2zHnQsiyNvlG5hTHMPRV2B3ALwG7yhs7Fh74GKuUS/uAG4ALgC8OLMesaYuY6sMY00OhX88j1lJLYzRCcxF1wSBxq+pdADXGt3bhARG5FzihyHQqIn8FvJ4Fibvr/XcsL3I13dDQBsC3nDFCTWG8nsXchA3xxtynAg+UPhcLD5y6/ru6vZY55S33YeqK0pWGyIUQE0xc36WHKO8S6RR3x8IDjT+r2aYt22tR1atVdUNVN7ZsyetBTi1wI386u+Wq+uoex30QOL30+TRWecwfXP9d3b5IhmZQDZUQokyfMNhY64qVLXNIK+1T3jk3wrFiy68HPiUiH2HlUDsTuEVVnxSR/xaRbcBXgV8D/jRSGWZDymmS2/J2x2CMc7UJe85irjJI3CLyBlbi3AJ8QURuU9XzOxYeuBz4S+CHWDnSFuNMi0EM6w3hJ7IUgiqXdaxGI4Vlf6bAkjUkRmG9fWeVTR2a6trNHVrOukbCZf86chC1JWuYEeUgiRwqX0hynLkVE7PcM6DPWHxMSx76nbfP8X3zqENejWab5bZkDTOgz6IHvt3X0ITyoLtMw3QhJ0G7YuKeCeXK6SP0psofauzbRV3DVPfbPhbYlSUKG6xbPlu2b9oZbRFCONryFsJz6e4W5aqWdQyWKOK2brmJOxPqRFVsD021UXERVXX/0OVaorDBvOWLoKlyx670rsePOVVyqcLuwiy3cYQuaxpTRL49AWOFecsNJ6YUlQk6PNYtN4xMMXEbRqaYuA0jU0zchpEpJm7DyBQTt2FkionbMDLFxG0YmWLiNoxMmU34qYgcAu6b6PQnAY9NdO4pWNr1wnyv+UdVtTYZ3WzEPSUisr8pfjdHlna9kOc1W7fcMDLFxG0YmWLiduPqqQswMku7Xsjwmm3MbRiZYpbbMDLFxG0YmWLibkBE3iciD4nIbet/ry19d6WIHBSRu0Xk/CnLGRoRuWB9XQdF5N1TlycWInKviNy+frb719ueKyJ7ROSe9f8nTl3OIdiYuwEReR/wP6r64cr2FwF/DZzDagXTm4CzSgsdzhYROQ74BrCd1XLLtwJvVNWvT1qwCIjIvcCGqj5W2vYh4DuqetW6YTtRVd81VRmHYpbbnx3Atar6fVX9FnCQldBz4BzgoKp+U1WfAK5ldb1LYQfwifXfnwBeP2FZBmPibucKEfk3Ebmm1EU7FXigtM+D6205kPO1VVHgSyJyQEQuXW97nqo+ArD+/+TJSheARWc/FZGbgOfXfPVe4GPA+1lVgvcDfwC8GZCa/XMZ2+R8bVVerqoPi8jJwB4R+fepCxSaRYtbVV/tsp+I/AXwd+uPDwKnl74+DXg4cNGmIudrOwpVfXj9/6Mi8nlWQ5Jvi8gpqvqIiJwCPDppIQdi3fIG1g+34A3AHeu/rwcuEpFnisgZwJnALWOXLxK3AmeKyBki8gzgIlbXmxUi8sMi8qzib+A1rJ7v9cDF690uBnZPU8IwLNpyd/AhETmbVbf0XuAyAFW9U0SuA74OHAbekoOnHEBVD4vIFcCNwHHANap658TFisHzgM+LCKw08ClV/XsRuRW4TkQuAe4HxlnBMBL2KswwMsW65YaRKSZuw8gUE7dhZIqJ2zAyxcRtGJli4jaMTDFxG0am/D+QVJX+7VBEUwAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
"# we load the subdomains from the hdf file\n",
"subdomains = MeshFunction(\"size_t\", mesh_brain_2d, mesh_brain_2d.geometry().dim())\n",
"hdf.read(subdomains, \"/subdomains\")\n",
"\n",
"# and plot\n",
"plot(subdomains)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Function Space"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "The [`FunctionSpace`](https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/functions/functionspace/FunctionSpace.html) class represents a finite element function space that defines basis functions (shape functions) over a discretized domain."
+ "The [`FunctionSpace`](https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/functions/functionspace/FunctionSpace.html) class represents a finite element function space that defines basis functions (shape functions) over a discretized domain."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"V = FunctionSpace(mesh_2d, \"Lagrange\", 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A function space is linked to a particular element type:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 12,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "\"FiniteElement('Lagrange', triangle, 1)\""
+ ]
+ },
+ "execution_count": 12,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"ele = V.element()\n",
"ele.signature()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The number of degrees of freedoms of each finite element depends on the element 'family' (e.g. Lagrange), degree of the element, as well as geometric dimensionality of the domain and dimensionality of the field of interest."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 13,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "dofs: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]\n",
+ "dof coordinates: \n",
+ " [[0. 0. ]\n",
+ " [0. 1.66666667]\n",
+ " [1.66666667 1.66666667]]\n"
+ ]
+ }
+ ],
"source": [
"cell = next(cells(mesh_2d)) # a cell from the mesh\n",
"print(\"dofs: \",V.dofmap().dofs()) \n",
" \n",
"print(\"dof coordinates: \\n\",ele.tabulate_dof_coordinates(cell)) # dof coordinates for that cell \n",
" \n",
"# change element type or degree"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A FunctionSpace over a mesh combines all finite elements of the domain:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 14,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[0. , 5. ],\n",
+ " [0. , 3.33333333],\n",
+ " [1.66666667, 5. ],\n",
+ " [0. , 1.66666667],\n",
+ " [1.66666667, 3.33333333],\n",
+ " [3.33333333, 5. ],\n",
+ " [0. , 0. ],\n",
+ " [1.66666667, 1.66666667],\n",
+ " [3.33333333, 3.33333333],\n",
+ " [5. , 5. ],\n",
+ " [1.66666667, 0. ],\n",
+ " [3.33333333, 1.66666667],\n",
+ " [5. , 3.33333333],\n",
+ " [3.33333333, 0. ],\n",
+ " [5. , 1.66666667],\n",
+ " [5. , 0. ]])"
+ ]
+ },
+ "execution_count": 14,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"V.tabulate_dof_coordinates()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function spaces exist:\n",
"- `FunctionSpace` for scalar fields\n",
"- `VectorFunctionSpace` for vector fields\n",
"- `TensorFunctionSpace` for tensor fields"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 15,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[0,\n",
+ " 1,\n",
+ " 2,\n",
+ " 3,\n",
+ " 4,\n",
+ " 5,\n",
+ " 6,\n",
+ " 7,\n",
+ " 8,\n",
+ " 9,\n",
+ " 10,\n",
+ " 11,\n",
+ " 12,\n",
+ " 13,\n",
+ " 14,\n",
+ " 15,\n",
+ " 16,\n",
+ " 17,\n",
+ " 18,\n",
+ " 19,\n",
+ " 20,\n",
+ " 21,\n",
+ " 22,\n",
+ " 23,\n",
+ " 24,\n",
+ " 25,\n",
+ " 26,\n",
+ " 27,\n",
+ " 28,\n",
+ " 29,\n",
+ " 30,\n",
+ " 31]"
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"U = VectorFunctionSpace(mesh_2d, \"Lagrange\", 1)\n",
"U.dofmap().dofs()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [`Function`](https://fenicsproject.org/olddocs/dolfin/1.3.0/python/programmers-reference/cpp/function/Function.html) class represents a function in a finite element function space."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"v = Function(V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Functions can be created from analytical [expressions](https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/functions/expression/Expression.html), by [projection](https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/fem/projection/project.html) into a function space."
+ "Functions can be created from analytical [expressions](https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/functions/expression/Expression.html), by [projection](https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/fem/projection/project.html) into a function space."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"scalar_field_expression = Expression(\"exp( - pow(x[0]-x0, 2) / pow(sigma_x, 2) - pow(x[1]-y0, 2) / pow(sigma_y, 2) )\", \n",
" x0=3, y0=3, sigma_x=1, sigma_y=1, \n",
" degree = 1)\n",
"scalar_field_discrete = project(scalar_field_expression, V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Their values are defined at the integration points and interpolated at all other points according to the basis functions (elements family, degree) of their function space."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 18,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Values at integration points: \n",
+ " [1.52299797e-08 2.08579001e-05 1.10431945e-04 2.26032941e-06\n",
+ " 2.08579001e-05 2.85655008e-02 1.51239760e-01 3.09558685e-03\n",
+ " 1.10431945e-04 1.51239760e-01 8.00737403e-01 1.63895538e-02\n",
+ " 2.26032941e-06 3.09558685e-03 1.63895538e-02 3.35462628e-04]\n",
+ "Value at arbitrary point (1.232, 3.434): 0.10287732297498928\n"
+ ]
+ }
+ ],
"source": [
"print(\"Values at integration points: \\n\",scalar_field_discrete.compute_vertex_values())\n",
"\n",
"print(\"Value at arbitrary point (1.232, 3.434): \", scalar_field_discrete(1.232, 3.434))"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 19,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 19,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbBklEQVR4nO2dX6hl113Hv79z5v7LJJORGGtsAhGUQCk00SEWBkRLrbEN1UcD7VPhvqikVCnGF1HwtdSHPnixQUtrSjENSJXqgA0hUNNmYtqmnSqlBowJDDGZzIy999xzz/n5cM5O9tl37b1/a++19vqzfx+4zJ179z1n3XPPZ6+1vvu31iZmhqIoeTIJ3QBFUfyhgitKxqjgipIxKriiZIwKrigZo4IrSsackRxERC8DuAFgAeCEmS/4bJSiKG4QCb7m15n5dW8tURTFOTpEV5SMIUklGxH9F4A3ATCAv2LmA8Mx+wD2AWCK6S/fQrc7bqoyFoio88+OpTLzBr/xOjPf2XacVPCfY+ZXiehnAFwC8AfM/Ezd8ecmd/D7z/ymVYMVBQAme7uhmxCc5eFR6zGX5k9clmRhojk4M7+6/vcqET0F4EEAtYIrSh0Sgcmj5CyQJzRtr5HkBFDQKjgRnQUwYeYb688/BODPxc+gZI3LHten2EM+hy+6nJwkPfi7ADy1nhedAfB3zPx162dSkiHEMPlt8XZ3Bn/utzmahXtuT7QKzsw/BvC+AdqieCTmuW0Ucsfw/HX0OPHYXAdXIiNmaaWckjvm3ymB+XsVFTwycpBWSlJyA8O3z8EJRQUfmDEJXMdG0FWRm3e2ArSoPzSbh26CERV8QFTumvl2Re7l7ubbcnJ0Mkjb+uDyxOTyZKGCD4TK3Twkr5O77mupYXOS4p0tZ5Kn/8olgMpdL3e55ytEXu7I35aTWfy9OyA/Sbkerajgnhm73JL5dvnNbyN3l+NjxOdJKv1XJ2Jyltu6IqxF7rKoi91p/waumR4tnD2WL5Y7Z7xJroJ7IkW5vZRxCsK0Qm6XYhf4eEyX+D4BqeAeiEnuoLXXgjDNJPdiZ5htCqaz5SDPExIV3DFDyJ3EggmLMK2QeyixC4Z+vjJDnVxUcId0lTsJYW2wCNOqci+2u2/24JLpcR4bR6jgjpDKnZ3MZSzm20C8cgP+2mJ14tjb7V2uqoI7QOVG5zCtKvdiJx7Jm5jO0ujhVfCBGIXcFmFaef6bmtxAt7aGOCmo4D0JvQVRcHqEaeVhcCFMSpLXEVPvroL3QOV2E6ZV5W6a/6YQflVPUmXhF9s06O+ggndktHJXdz3pEaaVRZDILfl+TMRwMlLBO5Cd3F22KupYvFKVe0Py4nvb9s0BgOlxt5/LGRXckiTk9r23WMfilTq5N4brHeXu+7MucX6i2d3pvC+bCm5BULlj2RDQcZhm6rWXnjZ1mcS56YpXVHAhXuSORVopPcK0pvl2Ibcvsd9uo+fHj/EEooILcCp3alIXOAjTmubbhXyLQC/PNL8t0QGo4K2MXu4BwrTQcrt4btsTxGJnMsiCExW8AZXbf5hWlXsZSVBWZpJwOq+C16By+yleMYVpMcsNNLcrdvlVcAMq93BhWlVu30GYK0IEarS3a30DQhW8gso9fJhWlXu5vaoAmxynU7UWKyp4iVHJ3fB7uAjTbOfbVbmrn8dC00lnuRXfpTIVfE0WcvcssgkRppWH5IXQyx2Z2JOZ9vBtqOCIXO6Byl5DhGnGXlsot+2xfUn1ZDJ6wYPIHbpWvYLrlWA2YZqp115u9Rd3Mk9DyMXu1OvWyaMWXOUetnhlKLldPg6QzsnCxGgFV7njCNMKuctC8rb/Ci86Drdl8pCMUnCne5cnKLfrMM1mvr36mlnuIcQukDxXrCeByd4uIEzrxYIT0RTA8wD+h5kf7ta08DjdATVhuV2HaTbz7dVzmuWmrfB3G+F5P7EX2/FsPmHTgz8K4AqAc57a4h2VO84wLSa5c0MkOBHdDeAjAP4CwKe8tsgTKnc8YZppvl3IPR1wmF5lEemQvA/SHvyzAD4N4La6A4hoH8A+AOzilv4tc4jK3X2Pcl9hWnkOHIPcofF1C+FWwYnoYQBXmfkyEf1a3XHMfADgAADOTe6IpsZwDHKXQ7M6Yg3TqnJvbfu5T3YT8+P2fm65w0kWu0h68IsAPkpEHwawC+AcEX2RmT/mt2n9SVVuibBSUgjTqnLvbJ+OiGfHiSwzi4xWwZn5MQCPAcC6B/8jlbsGwWO5lLcNX2Fa1+KVpvl2k9xNX+/DGE4aWV4HV7nTCdPKQ/JC4r0t+2H64TzLt3JvrF4VZn4awNNeWuIIlTudMM3Ua3eRu8vPDXFCGPo2RSayOu2p3O1yu14Jtvpas9yS+XZZ0L0tt1Uih3Pznks72/Psh+nZCD52uUNsq+QqTCvkdi12Qflx62TPlSwEj1XuoUg5TKvKfW7bzwbl148j32FnzXL3DCZH7i4VJi94zHIP0Xv7DtO6rgSzCdN8y+2S5TYntVdc0oKr3OG3VeobplXlvnXLreQ352n03L5IVnCVO58wzZfcSqKC5yp3OSgTHR9hmGYz3wZOy33bGbt9v+u4cRJPRrLYIUxnYS6XJSd46nLbSmx8jAB7lNuEabbzbddy27K1fSKqR0+RpH4rlbtbmOZzW6U+cpeH5IXc5xxIfj2i3js0yQiucucZppV77ULun9r6v5pX4B3enJ9tPSZleGcLNOtff5+E4Cr3cGGa6+KV1efNcpd7bYncdccNKX2MdzExEb3gTuWWEpHcIbZVWj3vMGGaSe7bp4fVl6GRtxZ7p75225mjqIK2UEQtuHO5Jb13onKnGKaZhuS2civNRCu4yp1GmGZbvFInd1ns89OfQMK1RVxbg7UR4nJZlIKr3Glsq9Rlvg3Uyy0VW5ETneAq97jCtKrc56ftIdu1hZ8wbbq9zG5n1agEH6vc5bk2MFyY5nolGNAut2m+bSP36NjbBQ671wZEI/hY5K7KXCbFbZX6hmlVue+YNEv+v8u8r3+7JgrBVW5/xSvl78UapknljpHlNjCJ5DZFJoILrnKPO0yryn1+Yl5Rdm3ZbdnnrVuzUS8ZDS64hDHJHese5a6KV0zz7arct5dyrrfGe7MTJwQVfPB7dCcid+zbKgFuwrTykNwkd64MudtqMMFTlltKH7lDVKYBdnum+ZL79skO3lqG2/yBtpa9byEcC3n8FhkSaoOAGKiTu0uCbjP/Lq6B28htE7AVf1Or3rvHJTJABU+aqaCTc33DPMnupJJVXRJZi/l3OWArilyKMtVioUnxnMVa8PJCk6LNxZbJxU0Pij3Ry5s9FHLTWvbJfPX6Fa9jseFieSVZ8XcovjY9Lv5dC204WU9n69HT0WL9+KuRUrGj6sZS0aPVE/Ba9qWF9Cp4ZEh6bskyRZudPyXVW5L9xPtutFAO1Kpyl08IUrnLvXeb3KbXoE3ucu/dJrep926VuyQyd+zJVfBIkAzbpo6vt0qGoq5v8dNWZvrWsl7u8s9K5S6PONrkrvbeQLvc5VGUVO6i9149fovcR5vDNJveG1DBg1P+Y9ceYzFnk/TuJOixbW7pI1l33bTy69py59R17jq5TY/TJnd59CGVu+i9gXa5N4brQrmL3htol7tr7w2o4MEo/4FdIAl7ym9aF0gCLNNmDHUUorfJXe29gXa5yyMRqdzl/KJN7vLoSip30XsDfuQGVHCvTGYnGx+2SObjroM2ye6iLoO24sP0PUAeqgHtcpdHJVK5y1lGm9ymkVab3Bu3KXIsN5BIJVsuFH/UpuvjwErs8nVw42PNNyvZ+rI4nmxcDzdxON/2doNAk+S+E/MCH4l5+eTcJzHvi/bgEeE6aJMk6a6DNkmSfm1xduNDgo/EvPy755CYm9AePAKms+XGQpO+SHp3Op5sLDTpy42T3dYbF1xb3GLctaVJ8nKoNrbE3AXagyeCpHd3HbRJknTboO3a4paNjyaqQ3MgrsS8IJbE3ERrD05EuwCeAbCzPv7vmflPnbZipEyPFhubPBiPEczHp7PNuvS+zI/PbGzyYOL68Y6T2/2aJD8//YmzxLzAR2K+MQSPIDE3IenBZwA+wMzvA3A/gIeI6P3OWzIiJIm661p0SZLuqqKtoCzlW4u9jY8mXCbmXcpQV5+v/rUpQw2dmJto7cGZmQHcXP93a/0x3pUQnpke88aa8CYkc+3JMW2sCTfB88nGmvC+XD/Zbb3HWJ3kxeozl4l5+cQVVWI+AKKQjYimAC4D+AUAn2Pm5wzH7APYB4BdpLVfdWpMjzd3dDHh+jLa4fzMxo4uJkxBm+l6eNPticrix5CYF3hJzC1DNdsyVUAYsjHzgpnvB3A3gAeJ6L2GYw6Y+QIzX9ii/G8ZszHccoCrklWb5YuuSlZNQVvT5bI352dPfZiOKT9OqMS8Sxnq6vGHT8xNWF0mY+ZrRPQ0gIcAvOSlRQlRlbzLtsmugjYbJnPa2HDRhCRoa8Mked3Q3SS5TWJe4CMx71KGCsjl9jX/BgQ9OBHdSUTn15/vAfgggB96a1HCTI5ONj4aj3UctLkqWbXZ+N9Ustq28OT6ye6pjyakiXmXMtTV19b/DpyYD4Wky7kLwN+u5+ETAF9h5q/5bVYeuOjhW5/DUdBmQ1vJap3kdYUwbZInn5gXDCw3IEvRvwvggQHakj3FH71NdEmS7jpokyTpTUFbIWF5b/QqJvGbqt+GTMwLvCTmHofgbWipakRISlZtLqNJkJSszo63Nm580IQpcLOVvu3xfCTmXcpQAX+JuStU8AiQBG02TI4372xiPMZR0Faeh9dVttWVszaJX/65PmWowHgScxMqeEAms5PWpaMFrkpWJzPa2ErZJabQramctam3d7lxA+A/MS+ISW5ABY+emNeGl6kL3fpI7zoxL/CRmDeVoYZEBQ/A5OhksKDt7ed0VLJaF7SZrknbSA9sih/TVsdSuYcuQ5Wgglco/5Fc3uVESgprw8vD5KbyVRvpgdPi+0jMu2zcsPo8bGLepUwVUMEbqZ6RfQpvE7RJendXQVtBXZJet9tLnfh1q9HqxI8pMX/7eSK8HFaHCm7BkMKbCL02vFqX3nTpzCS+bW9fEEtiLipDtcRnmSqggvfCNOeylV6SpLsO2iRJuiRoMy1EcSt9t8S8YPCtjiNEBXeMtJd3FbTZ0DVoM+1EWnd9vG71WZ34bRs6DrnVca/EPFJUcM/4Hta7Klm1DdpspAfse/vyzySZmEdwiQxQwQeHZvNWyV2VrEqCtibqVpbVDd3rbppg09sX0odMzAtirjGXooJHRCxrw9uuh5vEb5qvdx3iNyXmBcG3Oo4cFTwCXAVt7xzbv2TVJNOQ0idxc8AEGI/gh0fAXlxbSUmCNqvH67k2vJCnbi5edxeUOvH7DvGjvjlgIoxHcOD0HyeQ8JJ5eEGIteGmvdqaAjgfvX2UNwdMkHEJXiUS4U0MuTa8PNytq2yr26DRpre3ld5nYh7zVsdVupapAmMXvIrphRxYetdBm22SbioYaSpntentbYf4pnYFuzlgoqjgbQzUy9usDZcgDdoKmgI3n9IDzXc4jeLmgAmjgtsScFjvqmTVFLTV7bhaJ37dTQxthvhthTVRJOYe8V2HDqjg/ekpfIi14cBmONVUvmoS31Vv3zSvd5WYF3RKzDNABXeNcB7vqqLNhrrevSx7mTrxfUoPRLbVcRciKVMFVPBhsOzlfa4Nr24V3DScN4lvIz3QfYg/1sTcNSp4ojjbhNHwfnYlPdCvt48uMU+syAXIRfBiSLTrcKeDAPhcG266KWHd5TOT9EC9+D6kt0nMC8aemJvIQ/CC6twnEeFDrA0H7KQH7Hp7F/P6oW4OmDN5CV4lcuFdB22mJL3upoR1Q3ef0q8e37K393xzQCDf3hvIXfAqpnQzhPSOF76Yeve26+Em8W2kB/wO8U3PH2Vi7pk+ZarA2AQ3EWkv76Jk1bZ3tZEeGKC395yYjwEVvEpg4fsGbcUbv67oxbZ3DSF92+O5SMzHggrexkDDetdBm6maq6nSzaZ3dTGvBwT7tnu6OWAMDFGmCqjg3XDYy9usDW+i3GvVnQR8Sg+47e01MXeDCu4Cz8N625JV08YGNtKvjjd/fcghvo/E3DsRlakCKrgfJMN6QZIuDdoKmgI3G+lXx5/+2pBD/DEm5j6IWvDqPIUi2nHFmrL0LT1817Xhpq1/gXrxY5IeOC2+Jub9aX0XEdE9AL4A4GcBLAEcMPNf+m6YiayEN2CzCWNV5sbe2yC+jfSA3yF+o/SxJOYJ1qEDsh78BMAfMvMLRHQbgMtEdImZf+C5ba2YksjcpK9Sl6TbSNzpeI+9fdulMpvEXNmkVXBmfg3Aa+vPbxDRFQDvBhBccBMp9vK2JatV2Wp7V9she2RD/L5bHUvIuUwVsJyDE9G9AB4A8Jzhe/sA9gFgF7c4aJobohbeVdBmK5rH3t7VEN/0eNEl5gkgFpyIbgXwJIBPMvP16veZ+QDAAQCcm9xRv3ogMFEM649mnYI201C07vKZtWiRDfHLP2OTmBfkEK71rUMHhIIT0RZWcn+Jmb/a+1kjI+pevgUb6QHLa+QOhvhdpe+SmCunkaToBODzAK4w82f8Nyk8IYU3JenliqyCpmG7T+lXj28xZO8xr29KzAtSlHuoMlVA1oNfBPBxAN8johfXX/sTZv4nf82KiyGENwVtTdfDTdID9eLXJc02Q/wh5/U2ZahKPZIU/VkAbrcSwWp+MUloKFzG6Ty+JWgz9VBNRTA+e3sX8/rV48uOTy5Ui6xMFQhcyVYNEVIVHujQyzcEbUVAVFf0UvfGtuntYx/iJyd3pERVqmpKDVOV3sWw3pQEN1W62fT2Lob4vqWvPrfKbU9UgpvIpZcvC18ne7nooq7wxaf0gF1v72JeDzRdI68vQx2URMtUgQQEr5KL8G2YKqxspAfshvih5vVAvfh1j6XISU7wKlkIX/4dGtpvIz1g19u7mNcD7ob45Z/xlZjnXqYKZCB4leTn8XXDQcGwvsCV9EC4Ib5eDnNDdoKbiLaXt9kJRnhTQ6C+Z0phXg/otW7ATZkqMBLBq4QWng+PzEGb7QaPFtIDw8/rge6X7vrInUMduitGKXiVEMN68WW0uuKJOvE9Sg+EuXSXE0OWqQKAu5tPZ4arIZIUPjyy++PbVE0dHlld6qHZ3CqAmhydWPWak9lJp8tdLi+RiX8/6evmsIrN5XtPBU+ZCEsjfSGVW3KiCSW35ATuumNRwcdEogUbErmlo4gxyQ2o4EoDMVwnlsotYWxyAyq44hDb9LpNXpW7P5qiR0btJbRM6ROcqdztaA+eOrZBW6Lz8CoqtwwVXGkkhnl4FZVbjgquOMV3FZnKbYcKHiFDVzulQsxySwuVhi6gUsFzYAQFL7HLLWFouQEVfJxYvtFCz8NV7u6o4IpzXM7DVe5+qOCRovNwldsFKnguZDYPV7ndoJVsY0VwZ9M+DLHpgsrdjvbgiojQQVsVlVuGCq4kh8otRwWPGOugLbN5uAmV2w4VfMxE+IZsQuW2R0M2RYzpFsdDPrcIiWgO68qBeOUGtAdXEkDl7o724A1420bH4vLU2DaAqKJy90N78NzIKGhTufujPfjYsSx4ie16uMrdTGsPTkSPE9FVInppiAaNgZTeIFGjcrciGaL/DYCHPLdDaUAXnhhQuUW0Cs7MzwB4Y4C2KK7IaB5uROUW4yxkI6J9InqeiJ6fc5ovxmhJ6c2rclvhTHBmPmDmC8x8YYvGe1lH8YjKbY1eJgtE6m+cQZHeHVXlPoUKngijXXgSYFtjIA+5AdllsicAfBPAfUT0ChF9wn+zlMGJ8Q2tcvemtdCFmR8ZoiGKsoHK7QStZAvI8vDIqi7dO6m9uVXuVnQOnhBa8FJC5RahgudMLkFbFZVbjAqupIXKbYUKHpjc32BOUbmtUcETY7TzcJW7Eyp47uQwD1e5O6OCK3GjcvdCBVfiReXujQoeAbZvvFHMw1VuJ2gl2xjIYR5uQOVuR3twJUlUbhkquJIcKrccFTwS9M0oQ+W2Q+fgiTLGO55I5FaxN9EeXEkClbsbKrgSPSp3d1TwiNA36WlU7n6o4AmTe8GLyt0fDdkSJ3fJm1C529EeXEkSlVuGCq4kh8otRwWPDH3zNqOvjx0quJIMKrc9KriSBCp3N1RwJXpU7u6o4BGib+h30NeiHyq4Ei0qd39UcCVKVG43aCVbpOgbXHGB9uCKkjEquKJkjAquKBmjgitKxqjgipIxIsGJ6CEi+g8i+hER/bHvRimK4oZWwYloCuBzAH4LwHsAPEJE7/HdMEVR+iPpwR8E8CNm/jEzHwP4MoDf9tssRVFcICl0eTeA/y79/xUAv1I9iIj2Aeyv/zu7NH/ipf7NG4SfBvB66EZYkFJ7U2orkFZ775McJBGcDF/jU19gPgBwAABE9DwzX5A0IDQptRVIq70ptRVIq71E9LzkOMkQ/RUA95T+fzeAV7s0SlGUYZEI/m0Av0hEP09E2wB+F8A/+G2WoiguaB2iM/MJEf0+gH8GMAXwODN/v+XHDlw0biBSaiuQVntTaiuQVntFbSXmU9NpRVEyQSvZFCVjVHBFyRingqdU0kpEjxPRVSKK/no9Ed1DRN8goitE9H0iejR0m5ogol0i+hYRfWfd3j8L3aY2iGhKRP9ORF8L3ZY2iOhlIvoeEb3YdrnM2Rx8XdL6nwB+A6tLa98G8Agz/8DJEziGiH4VwE0AX2Dm94ZuTxNEdBeAu5j5BSK6DcBlAL8T8WtLAM4y800i2gLwLIBHmfnfAjetFiL6FIALAM4x88Oh29MEEb0M4AIztxbluOzBkyppZeZnALwRuh0SmPk1Zn5h/fkNAFewqjCMEl5xc/3frfVHtGkuEd0N4CMA/jp0W1zjUnBTSWu0b8JUIaJ7ATwA4LmwLWlmPeR9EcBVAJeYOeb2fhbApwEsQzdECAP4FyK6vC4Rr8Wl4KKSVqU7RHQrgCcBfJKZr4duTxPMvGDm+7GqfHyQiKKcBhHRwwCuMvPl0G2x4CIz/xJWKzx/bz3dNOJScC1p9ch6LvskgC8x81dDt0cKM18D8DSAhwI3pY6LAD66ntd+GcAHiOiLYZvUDDO/uv73KoCnsJoeG3EpuJa0emIdWn0ewBVm/kzo9rRBRHcS0fn153sAPgjgh2FbZYaZH2Pmu5n5Xqzes//KzB8L3KxaiOjsOmgFEZ0F8CEAtVeCnAnOzCcAipLWKwC+IihpDQYRPQHgmwDuI6JXiOgTodvUwEUAH8eqd3lx/fHh0I1q4C4A3yCi72J14r/EzNFffkqEdwF4loi+A+BbAP6Rmb9ed7CWqipKxmglm6JkjAquKBmjgitKxqjgipIxKriiZIwKrigZo4IrSsb8P+nGjzaM/WqbAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
"plot(scalar_field_discrete)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similar for Vector fields and tensor fields"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"vector_field_expression = Expression( (\"cos(x[0])\", \"sin(x[0])\"), degree = 1)\n",
"vector_field_discrete = project(vector_field_expression, U)"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 21,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 21,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVnUlEQVR4nO3de3SV9Z3v8fc395AQEnKBSICAXLQgBQ2CjTLnUOVYbR1rp/Zmtc50WOOZrtrV6ellOl3TnjlrndOzOl2d03VOp9jaVq162qqnFbEttjCIIki4yJ0gcofsJFxyIeT6PX8kRmCg2dnszbP98XmtlZU8ez95nk9gf/bvueWJuTsiEo6MqAOISHKp1CKBUalFAqNSiwRGpRYJTFYqFlpWVubV1dWpWHQw2rq6eOvkCaaXlpGTmRl1HHmXqaura3L38gs9l5JSV1dXs379+lQsOhg/rFvHt195mU/dMJev1C6IOo68y5jZ/os9p83viGyNNQDwq+1b6ezpiTiNhESljsjWWAyA5o4OfvdmfcRpJCQqdQRaOjvZf+rk4PSTW96IMI2ERqWOwLZYA2UjRpCXlUVVUREZZtQ3N0cdSwKhUkegdMQIlt/3IOOLRlGUk8uTH7mXySUlUceSQMR19NvM9gGtQC/Q4+41yQ7S2tFJYV4OZpbsRaedaaVlAJQXFAyO0JkZen+V5BjOK+k/uvvsVBQaYPu+Y3zoG4/yP//vCtZs309Xd/hHhMtHFNDccZrevr6oo1xWB2Inh55JEpaS89TxWLFpD6/vOjg47cDxltM8vXITT6/cxIjcbG56z0QWXDeZm2dOomTkiKiipkxFQQF97hzv6KC8oCDqOCn31rHj/OD5VzEzvv3ZO2k9fYadBxvZdTDGnfPfQ0lhftQRgxBvqR34vZk58EN3X3L+DGa2GFgMMGHChCEXuOnNIzy9ctNFnz/d2c0fNu5h96Em9hxp5v7bbqC0KKwX/sdnzmLR1VMoys2NOkpKHWluYckLr7H0te30uTOhopi7vvEoh5pODc5TPXY0N8+cFGHK6C1/8Q2e+MnLfP+RBykalfggFm+pa939iJlVAMvNbKe7rzp7hoGiLwGoqakZ8s4LD33ofXz2jnmD0yfbOvjYPz1OV08vs6dcxYLrJrPguslUjx0d/0/zLjOpuIRJxeEeIGtuaefR367jVy9vobund/Dxw02nmDqunA/XzuSa8RVcM6GCqeMueMXjFSUzM4Mjh06wfu1eFi6amfBy4iq1ux8Z+Bwzs+eAG4FVf/q7/rS8nCzyzlr9zgMx/uG+W6mdMYlRBXmXsmhJA0eaW3hk2WsciJ1kTEkhTSfbOTNwnCQzI4Nv3Hcr104YE3HK9FIz/2oyMozXXqlPbanNrADIcPfWga8XAf814TVexNzp45O9SInQVaVF/OOnFw1OuzttZ7poPNlG06l2Tnd2R5guPRUV5TPjuipef+1Nurp62LntMLPmTBz2cuIZqccAzw2casoCnnT33w57TXJFMzNG5ucyMj+XyZWlUcdJO319zt49DcyZO4ktmw/yxYceY9o1lakptbvvBd6bSFARiU9GhrHsNxt5/tk6AHbtOMK0aysTW1Yyg4lI4j770EIqxhQNTid6IZZKLZImRhTk8oWv3Dk4nejFlSq1SBqZO/9qFt0xCwBDI7VIEP7m87cxurSQBDutUoukm5FF+Xz+Sx9IeJ86smu/ReTiav9sOsWjE7tUVCO1SJqacV1iF2Sp1CKBUalFAqNSiwRGpRYJjEotEhiVWiQwKrVIYFRqkcCo1CKBUalFAqNSiwRGpRYJjEotEhiVWiQwKrVIYFRqkcCo1CKBUalFAqNSiwRGpRYJjEotEhiVWiQwKrVIYFRqkcDEXWozyzSzjWa2NJWBROTSDGekfhjYkaogIpIccZXazKqAO4EfpTaOiFyqeEfq7wFfBvouNoOZLTaz9Wa2vrGxMSnhRGT4hiy1mX0QiLl73Z+az92XuHuNu9eUl5cnLaCIDE88I3UtcJeZ7QOeBhaa2RMpTSUiCRuy1O7+NXevcvdq4OPAH939vpQnE5GE6Dy1SGCyhjOzu68EVqYkiYgkhUZqkcCo1CKBUalFAqNSiwRGpRYJjEotEhiVWiQwKrVIYFTqCL2xcX/UESRAKnVEThxv46dLVkYdQwKkUkekftcxtmw+yFtvxqKOIoFRqSNSv/MoAC/8ekPESSQ0KnVE6ncfA2D5i1voON0VcRoJiUodkbdH6tPtnax4aVvEaSQkKnUETp08TayhZXB66XN1uHuEiSQkKnUE6ncdZUJ1Gfkjchg/oZQJE8vYt1c3a5TkUKkjUDWhlH/92V8ztrKY7JxMvvrNu5l0dUXUsSQQaVPqho6TdPX1RB3jshhbWUx2diajSws53tQWdRwJTNqUurGzlUV/+G98eeMTPH+ojubO1qgjpdzo0gJOnjxNT09v1FEuG3dnw9EjUccI2rDuUZZMSw9vYFVs+zmP9XofKxu2s7JhO4YxY1QVN1dcwy0V1zClcCxmFlHa1BhdOhKAkyfaKSsvijhN6q05eIDvrFnNxFHFzBozlr0njrO9McbWWIy/vr6GMYWFUUeMXHdvL9mZmZe0jMhKfbC9iXVNe855rKvvnRHLcXa0HCY3M4vcjGxKc0YyOjes//SPfPxG7rrnBkpGh/VznW9zwzH++dXVrD7Yf637wZZTzPrX73Om553drZuqxl/xpX5s80a+8+pqXrr/QSoKEv+3iKzUD01bxEPTFg1ON5w5xT3/9h3ys3K4qWwat1Rcw01l0xiZnR9VxJQLvcz1zc1897VX+N2b9ec83trZxS0TJjKjooKZ5WOYUVHBmEt4EYeioqCQtu4uVux7i4/NuC7h5URW6vPFzpzi+3MfZFbxRLIyLm3zQ6J3uKWFZ3ZuIyczk3njqoi1txNrb6O9u5vuvl4emD2H2vETo46ZVm6eMJHsjAxW7NsbRqmvK54QdQRJonFFRXy1dsG/e7ytq4tYextdvVfOwcF4FebkMG/ceFYf2E9LZydrDx3ktqunDHs5aXP0W64MhTk5TC4ZzTVl+iOKZ+tz59/2vcV7Kio43d3NXU89zqoD+xJalkotkgYyzKg7eoQlda8DcKDlFIme61GpRdLE386dx9TRpYPTiZ7CValF0kRuVhbfvvU/kTFQZo3UIgGYPbaSB2dfD2ikFgnGF+fXMnFUcepGajPLM7N1ZrbZzLaZ2bcSXJeIxCE/O5v/8f5FCY/U8Zyn7gQWunubmWUDq83sRXd/LaE1isiQ5lWNpzg/sasphyy199+S4+3fD8we+NBtOkRSbHppWULfF9c+tZllmtkmIAYsd/e1F5hnsZmtN7P1jY26i4dIVOIqtbv3uvtsoAq40cxmXmCeJe5e4+415eW6WkgkKsM6+u3uJ4GVwO0pSSMilyyeo9/lZlY88HU+cCuwM9XBRCQx8Rz9rgR+ZmaZ9L8J/MLdl6Y2logkKp6j328Acy5DFhFJAl1RJhIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosEZshSm9l4M1thZjvMbJuZPXw5golIYob8o/NAD/B37r7BzEYCdWa23N23pzibiCRgyJHa3Y+6+4aBr1uBHcC4VAcTkcQMa5/azKqBOcDaCzy32MzWm9n6xsbG5KQTkWGLu9RmVgg8A3zB3VvOf97dl7h7jbvXlJeXJzNjkNydl/buiTqGBCiuUptZNv2F/rm7P5vaSFeGY21t/MvaNbh71FEkMPEc/Tbgx8AOd/9u6iNdGbY1NrCtMcYbDceijiKBiWekrgU+DSw0s00DH3ekOFfwtsQaAHhiy+aIk0ho4jn6vdrdzd1nufvsgY9llyNcyLbGYgAs3b2LU2fORJxGQqIryiKyrbF/pO7s7eGZHdsiTiMhUakjEGtvI9bePjj95NbNOmAmSaNSR2BLrIH548ZTlJvL9NIybps8hd3Hm6OOJYFQqSPw3jGV/Pyej1I1sgh35yu1C5heWhZ1LAlE2pR679Fmmk61Dz1jAMpGjMDMKC8oIHb6yviZ5fKJ5xc6Lgt35/avPcK1EypYMGsyC66bzLSqcvpPk4epvKCAk2fO0NnTQ25W2vxXpJR7F3S+jOW9P+oowYrslfTUHzeybN2Ocx7LzMxg2/4Gtu1v4AfPr2FMSSG3zJzMglmTmTt9PLnZYb3wy0cUAHC8o4PKkSMjTpNa7r1wZine9r8gZz7kzIHu7dC9De/ZjhX+FyyrKuqYkTt5op3ikoJLWkZkLenq6aW9s/ucx84/Anymq4fTnV20n+miq6c3uFJ/bu58vjDvfWRnZkYdJWXcHTpfwtu+Bz31/Q92xPCOX541Vzbk/wVc4aV++rFXePSHK3jsV59jbGVxwsuJrCUPLKrhgUU1g9P7G07wkW/9jEljR7Ng1mRumTmJWZOvIiszbXb7ky4/OzvqCCnlna/ibd+F7jfOfcIyIf+TWNYMyJ4BWVMwy4kmZBqZPHUM7rD21T38+Udqhv6Gi0iboa+zu4fnvvUZxpcn/g4l6cN7D0PPbsi5ETKroS8GvTHoawRvxbJvxPJ1tfHZZl9fTV5eNmtfrQ+j1NOq9OuaIbHMcVDwGS50mNP7ToO3XvZM6S4nN4vZNdXUrdtLY6yFjevfYtEd7x32csLdtpW0ZRkjsMwxUcdIK319zu+Xbaa8oojurl7+5oFH2LntSELLUqlF0kBGhtHacobnn60DoOVUBxfczIlnWUnMJSKX4O6PzuXaGe/c/i/RazRUapE0kZmZwd99/UNkZ/ef4kz0uiuVWiSNTKwu476/vAUAS3D7W6UWSTP3fuomrp46RvvUIqHIysrkS3//zmb4sL8/yXlEJAmmTB9LYVFeQt+rkVokTSV6/bdKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlVokMPH80flHzSxmZlsvRyARuTTxjNQ/BW5PcQ4RSZJ4/uj8KuD4ZcgiIkmQtH1qM1tsZuvNbH1jY2OyFisiw5S0Urv7Enevcfea8nLdw1skKjr6LRIYlVokMPGc0noKWANMN7NDZvZXqY8lIoka8h5l7v6JyxFERJJDm98igVGpRQKjUosERqUWCYxKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlVokMCq1SGBUapHAqNQigVGpRQKjUosERqUWCYxKLRIYlTpC7h51BAmQSh0Rd+elx1dFHUMCpFJH5Ni+GE/992c1WkvSqdQRqa/by8FdR9i0YmvUUSQwKnVE6uv2ArD0h8sjTiKhUakjUr/xLQBeeW4dzUdPRJxGQqJSR8DdB0fq3p5efvvoHyNOJCFRqSMQO9BES3Pr4PSyR16it7c3wkQSEpU6AvUb9vL+T93CqLKRTL1hMh/90l0c2nUk6lgSiCH/lK0k3w2L3svNH57Hf675MmfaO7n7cx+IOpIEJK6R2sxuN7NdZrbHzL6aiiC7dhxh986jV8R52/yCPABGV5Zw4tjJiNNIaIYcqc0sE/jfwG3AIeB1M/uNu29PZpBRxSP4zL3/h+KSAubXTmH+zdOYU1NNbm52MleTVkrGFNN2sp3Ojk5y83OjjnNZtLZ08MqqXdz+wdlRRwlWPJvfNwJ73H0vgJk9Dfw5cEmlfuInL/P8c3XnPOY4zU2tvPDrjbzw643k5mYxp2YS82uncvN/uIZRxSMuZZVpp7yqlNFji2k90R58qTs6uvh/v3ydX/x8DX+28FqmTB1L/e6j7NndwJ5dx/jS1z/E+ImlUceMlLtTf7iJqePKMLOElxNPqccBB8+aPgTMO38mM1sMLAaYMGHCkAstKSmgelL5OY9tbTlDV2/P4HTlVSVUTy6n+upyCkfmxRH13eX+b97L/d+8N+oYKdXd3csLv97Akz9dzYnj7QCDb9pvKysfyfHm1iu+1D94fg0/enEtz/7jA1SPHZ3wcuIp9YXeMv7djq+7LwGWANTU1Ay5Y3zn3ddz593XD07v2XWMzy/+CdfP7R+Z59dOpXJcSRzxJB319vbxh99t4fEfr+LY0VPnPJeXn80n7q9l6vRKpkwbS8nogohSppe508fzoxfX8vLWt1Je6kPA+LOmq4Ckn38ZUZDDL5d9kYKCsDdDrxTHm9vIy8vhno/No7mpjePNrf2fm9pobmqlsDCPufOvjjpmWpk95SoK83J4ectePn3rDbh7Qpvh8ZT6dWCqmU0CDgMfBz457DUN4aqqxN+ZJP2UVxRRvrDoos93d+tim/NlZ2Zy04xqVmzcQ/3hJl7fdYBPLrx+6G88z5CntNy9B/gc8DtgB/ALd9827DWJnCU7OzPqCGmlr895/KU6crIy6enr4/5vP8n+hsROd8Z18Ym7LwOWJbQGERlSRoZRUpjPC2t3ANDZ3UuiB8B1mahImrhz3rW87z3Vg9OJntZSqUXShJnx9U+9n4K8nP7pBJejUoukkcrRRTz84ZsBtPktEop7bp5FzbQqEh2rLRW/QGFmjcD+IWYrA5qSvvLkUb7EpXM2SO988Wab6O7lF3oiJaWOh5mtd/eaSFYeB+VLXDpng/TOl4xs2vwWCYxKLRKYKEu9JMJ1x0P5EpfO2SC9811ytsj2qUUkNbT5LRIYlVokMJGU+nLcyDBRZvaomcXMLO3+yJWZjTezFWa2w8y2mdnDUWc6m5nlmdk6M9s8kO9bUWc6n5llmtlGM1sadZbzmdk+M9tiZpvMbH3Cy7nc+9QDNzLczVk3MgQ+kewbGSbKzBYAbcBj7j4z6jxnM7NKoNLdN5jZSKAOuDuN/u0MKHD3NjPLBlYDD7v7axFHG2RmXwRqgCJ3/2DUec5mZvuAGne/pAtjohipB29k6O5dwNs3MkwL7r4KOB51jgtx96PuvmHg61b6f799XLSp3uH92gYmswc+0uZIrJlVAXcCP4o6SypFUeoL3cgwbV6Y7xZmVg3MAdZGm+RcA5u3m4AYsNzd0ynf94AvA31RB7kIB35vZnUDN/JMSBSljutGhnJxZlYIPAN8wd1bos5zNnfvdffZ9N/L7kYzS4tdGDP7IBBz97ohZ45OrbtfD3wA+NuBXcFhi6LUl+VGhqEa2Fd9Bvi5uz8bdZ6LcfeTwErg9oijvK0WuGtgv/VpYKGZPRFtpHO5+5GBzzHgOfp3VYctilIP3sjQzHLov5HhbyLI8a4zcCDqx8AOd/9u1HnOZ2blZlY88HU+cCuwM9pU/dz9a+5e5e7V9L/m/uju90Uca5CZFQwc/MTMCoBFQEJnYC57qdP9RoZm9hSwBphuZofM7K+iznSWWuDT9I8ymwY+7og61FkqgRVm9gb9b97L3T3tTh2lqTHAajPbDKwDXnD33yayIF0mKhIYXVEmEhiVWiQwKrVIYFRqkcCo1CKBUalFAqNSiwTm/wPffWKEojUgswAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
"plot(vector_field_discrete, wireframe=True)\n",
"# for this to look prettier:\n",
"# - define a domain centered around (0,0)\n",
"# - increase the number of elements"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving PDEs "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Above, we have explored the basic components that are needed to define and discretize functions over arbitrarily shaped domains.\n",
"We are now going to see how to use these tools to solve partial differential equations with given boundary and initial conditions over these domains."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a diffusion problem, such as the transport of heat in an isotropic homogeneous medium:\n",
"\n",
"$$\\frac{\\partial \\Theta}{\\partial t} = - \\alpha\\, \\nabla^2 \\Theta\\, ,$$\n",
"where $\\Theta$ is the temperature and $\\alpha$ a real coefficient called (thermal) diffusivity.\n",
"\n",
"With \n",
"- initial condition $\\Theta(x, t=0) = \\Theta_0$\n",
"- flux boundary condition $\\frac{\\partial \\Theta}{\\partial n}$ given"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To solve this equation, we start by writing the weak form with test function $\\nu$:\n",
" \n",
"\n",
"$$\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right )\\, \\nu\\, dV= \\alpha\\, \\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dV \\tag{1}$$ \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the chain rule theorem, the second term of above equation can be written as:\n",
"\n",
"\\begin{align}\n",
" \\alpha\\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dV \n",
"&= \\alpha \\int_{\\Omega} \\nabla \\cdot \\left( \\nu\\, \\nabla \\Theta \\right)\\, dV - \\alpha \\int_\\Omega \\left(\\nabla \\nu \\right) \\cdot \\left( \\nabla \\Theta \\right)\\, dV \n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Applying the divergence theorem to the first term yields:\n",
"\n",
"\\begin{align}\n",
"\\alpha \\int_{\\Omega} \\nabla \\cdot \\left( \\nu\\, \\nabla \\Theta \\right)\\, dV\n",
"&= \\alpha \\int_{\\partial\\Omega} \\nu \\, \\left( \\nabla\\, \\Theta \\right) \\cdot \\mathbf{n} \\, dS \n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Equation (1) now reads:\n",
"\n",
"\\begin{align}\n",
"\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right )\\, \\nu\\, dV= \\alpha \\int_{\\partial\\Omega} \\nu \\, \\left( \\nabla\\, \\Theta \\right) \\cdot \\mathbf{n} \\, dS - \\alpha \\int_\\Omega \\left(\\nabla \\nu \\right) \\cdot \\left( \\nabla \\Theta \\right)\\, dV \\tag{2}\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\mathbf{n}$ is the outward facing unit vector on the domain boundary $\\partial\\Omega$. \n",
"By construction $\\nu \\equiv 0$ on the parts of the boundary on which essential boundary conditions apply.\n",
"\n",
"Therefore, the boundary integral only yields non-zero values where natural boundary conditions are applied.\n",
"These correspond to the von-Neuman boundary conditions that impose a value on the heatflux $q = \\alpha\\left(\\nabla \\Theta \\right)\\cdot\\mathbf{n}$. \n",
"\n",
"Separating terms that depend on $\\left(\\Theta,\\nu\\right)$ and $\\left(\\nu\\right)$, respectively, yields:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right )\\, \\nu\\, dV + \\alpha \\int_\\Omega \\left(\\nabla \\nu \\right) \\cdot \\left( \\nabla \\Theta \\right)\\, dV= \\int_{\\Gamma_{vN}} \\nu \\, q \\, dS \\tag{3}\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We approximate the time derivative by a backward difference scheme:\n",
"\n",
"$$\\left(\\frac{\\partial \\Theta}{\\partial t}\\right)^{k+1} \\approx \\frac{\\Theta^{k+1} - \\Theta^{k}}{\\Delta t}$$\n",
"\n",
"Equation (3) then becomes:\n",
"\n",
"\\begin{align}\n",
"\\int_\\Omega \\frac{\\Theta^{k+1} - \\Theta^{k}}{\\Delta t}\\, \\nu\\, dV + \\alpha \\int_\\Omega \\left(\\nabla \\nu \\right) \\cdot \\left( \\nabla \\Theta^{k+1} \\right)\\, dV= \\int_{\\Gamma_{vN}} \\nu \\, q^{k+1} \\, dS\n",
"\\end{align}\n",
"\n",
"Including the initial condition $\\Theta_0$ and separating terms that depend on the unkown $\\Theta^{k+1}$ yields:\n",
"\n",
"\\begin{align}\n",
"\\Theta^0 &= \\Theta(t=0) = \\Theta_0 \\tag{4}\\\\\n",
"\\int_\\Omega \\Theta^{k+1}\\, \\nu\\, dV + \\Delta t \\, \\alpha \\int_\\Omega \\left(\\nabla \\nu \\right) \\cdot \\left( \\nabla \\Theta^{k+1} \\right)\\, dV &= \\Delta t\\, \\int_{\\Gamma_{vN}} \\nu \\, q^{k+1} \\, dS + \\int_\\Omega \\Theta^{k} \\, \\nu \\, dV \\tag{5} \\\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To translate this into FEniCS, we first need to define a mesh, a function space over this mesh and suitable functions:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# mesh & functionspace\n",
"mesh = RectangleMesh(Point(-5,-5),Point(5,5), nx=10, ny=10)\n",
"V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
"\n",
"# functions\n",
"T = TrialFunction(V) # temperature (k+1)\n",
"v = TestFunction(V) # test function\n",
"T_k= Function(V) # temperature at previous time point (k)\n",
"\n",
"# other parameters\n",
"q = Constant(0) # assuming isolated boundaries, zero heatflux in/out of domain\n",
"alpha = Constant(1) # arbitrary value for thermal diffusivity\n",
"delta_t = 0.2\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This allows us to encode equation (5):"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
- "LHS = T * v * dx + delta_t * alpha * inner(grad(v), grad(T)) * dx\n",
- "RHS = delta_t * v * q * ds + T_k * v * dx"
+ "a = T * v * dx + delta_t * alpha * inner(grad(v), grad(T)) * dx # a = a(T, v)\n",
+ "L = delta_t * v * q * ds + T_k * v * dx # L = L(v)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additionally we need to specify initial and boundary conditions to find a unique solution to our problem."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# initial value\n",
"T_0_expression = Expression(\"exp( - pow(x[0]-x0, 2) / pow(sigma_x, 2) - pow(x[1]-y0, 2) / pow(sigma_y, 2) )\", \n",
" x0=0, y0=0, sigma_x=5, sigma_y=5, \n",
" degree = 1)\n",
"T_0 = project(T_0_expression, V)\n",
"\n",
"# Natural Boundary Conditions are encoded directly in the variational form\n",
"# -> heatflux q\n",
"\n",
"# Essential boundary conditions can be specified separately:\n",
"# def boundary(x, on_boundary):\n",
"# \"\"\"\n",
"# identifies boundary of the domain\n",
"# \"\"\"\n",
"# return on_boundary\n",
"# bc = DirichletBC(V, u_D, boundary)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now solve the time-dependent PDE:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 25,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Time step 00, t=0.00: T(0,0)=0.97\n",
+ "Time step 01, t=0.20: T(0,0)=0.94\n",
+ "Time step 02, t=0.40: T(0,0)=0.91\n",
+ "Time step 03, t=0.60: T(0,0)=0.89\n",
+ "Time step 04, t=0.80: T(0,0)=0.86\n",
+ "Time step 05, t=1.00: T(0,0)=0.84\n",
+ "Time step 06, t=1.20: T(0,0)=0.82\n",
+ "Time step 07, t=1.40: T(0,0)=0.80\n",
+ "Time step 08, t=1.60: T(0,0)=0.78\n",
+ "Time step 09, t=1.80: T(0,0)=0.76\n"
+ ]
+ }
+ ],
"source": [
"# Create VTK file for saving solution at each time step\n",
"vtkfile = File('output/temp_diffusion.pvd')\n",
"\n",
"# Compute solution\n",
"T = Function(V)\n",
"T.rename(\"T\", \"Temperature\")\n",
"t = 0\n",
"T_k.assign(T_0) # assign initial condition\n",
"for k in range(10):\n",
" \n",
" # increment time\n",
" t = t + delta_t\n",
" \n",
- " # solve variation problem\n",
- " solve(LHS==RHS, T)\n",
+ " # solve linear variation problem\n",
+ " solve(a==L, T)\n",
" \n",
" # save to file\n",
" vtkfile << (T, t)\n",
" print(\"Time step %02d, t=%.2f: T(0,0)=%.2f\"%(k, k*delta_t, T(0,0)))\n",
- " # update T_k\n",
+ " # update T_k with current solution\n",
" T_k.assign(T) "
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 26,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Temperature at center at t=0: 0.9999999999999999\n",
+ "Temperature at center at t=1.80: 0.7612013861310724\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD4CAYAAAAn+OBPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2db8ht2V3fv79z7vPnOjN3JnTSBjJDJ/QPNKhgGUJpXlhMlBiH+KYULYrgi6FQIQFDapzXfSWoBYVysZQWAyKoWKQSJ7S+6IuINzGxpKMSpOrEiJpq7mRynz/3Oasvzt3nWXvt9ef3+63f2nud5zlfeLj3nLPW2uvssz/r92etvTY553DQQQfdHq2W7sBBBx00rw7QH3TQLdMB+oMOumU6QH/QQbdMB+gPOuiW6c4SBz2mU3d39TSwMhxzVpR4P3IMipd1sTYib03KBS9drP1VWCbahfj7idOUaiPW55KSbQlE2omgRL1kext++cl7m/DzSKXgLdqUy6TLRd7bRL5ArG5MsboRPdx87a+dc++MfbYI9HdXT+Ofv/Nf2TZ6epI42OnkLXdyNHlvcxo/FZuT8ftXp+tJmauTMZVXxxR8PiUq+t5xnLyr40i/pl8BV4lTUNIm0n6NVhfyOuvzRFuXkbKJ9tcXU3DW5y77OlZvfT4Fa312Ne7X+eNJmdXZ9D06j3yBR2fT9wDgLHESArlUfU+ffvif/yT12TLuvaWFB0TASxQCH1MIPEcx4JNlmcBrtDm2B75lu4Ni52T7fvm8cs495zeNXRsxwxEzMMnrMnUdB6LK63r/Y3rmiRoU/RGYiln5SRnGhcetl7q4o/WFVr4llJpjpPqfGuC454YFufI3a6IZwN9v6HMnSHBSUq59SRyLEF50Eisfk4WVnwN4/1jc42nDk1EbRtZ+Uicy4HOtvViNwd9f6IUWHpBZeY5rH6p3K691uzdH9YNNjcs/h7Wf1rFDQ+TiD1Jc31ztJ/SlE9LAynNc+0kdY9eyBjwL2K3gz0nq5kfbUAy+Xbn4AtHdU7HF3z/oG46Ag5ZO4Fla+RrrnvusdgCyCjG01l6T0Ktx8Vtbewn4+wV9BfA1CTyOWlsXDWStXflWVv9g7T01AH9/oOd++UrXPjZyW7j2c1r5FtY9V2fusCOU5JyN6s2c0FMbHmPw9wP6Spfe2spbJHlqrXwOeI2qrbZBvO/LwtpP21w2oRcV10gZgt8/9BLgF7LynBV4HGkt1qClMvO17Und/GjZSBs3xsU3Vt/Qz5C0k8oigdfKyktlCfu+tG2xbmL2hN6g0xMTJvqF3gh4yTp7jawsw5xW3tq6WxynJ2s/7cMMmEim3SrZ6BN6zZdqsM5+zgRetE8NrPwcsGuPKR24rKTxxJok9GYCvz/oDYHflwRe7U01XFh0035u9KdVjXdRa+2n7XUct88Afl/Qz2Dhe0zgcRW7+DnA65Jrccgt4M8f166tQa1cfK61F0sKvpCbfqBfwKWXqFUCr+Wts9q2OFDPbfWtrf20fX2ytSTVnL302hbw0wf0DYDnJvCsrHxrSa281p2WwNzK6s9h7bt28Qc1At8MeiJaE9HvEdFviCp2buG5KiXw5rTyNbG76ng1dSVLaxtb+2kb4W+qd/HVK/QagG9p6T8K4A1RjUbAW1v5FlM2tcttU9awtXWfo51tW4KyRuFQt9bfGHyTq5mIXgDwfQB+gV1pjy18yVXUXjy1U3RWyboaadpsZe0ndQ1c/FkTer4Mr30rE/azAD4BpPYpBYjoVSJ6QEQPLlx5Y7+R7p6yv3RrK6+5cCZ1Kje6jFlB/lx4/fSb9DisxGCk/6kbh6KeTySHwVmsowrLlOCnXHz21PLAwdJ75BHRKwD+0jn3uVw559x959zLzrmXj9d3eY0Lv+DcwMdUumhqgI9BwEnYaUAfjpX6k4ozCKS+Swr8FPyjch2CXw0/UAW/hQ/yfgAfIaIPAzgFcI+IftE590Oq1hRfJHWyWgMvvVjCCzC5DRQzhk9nv7UJOX1Z6bbXfh9XF+PztDmabn09HC88ztXJdPvssP7V8XTb7KtjGm19fXVCo+2xJ6+D8tsyq8l22Ven68l22ZuTO5Mtszend6JbZvvXcnT77FA+L4ytsQEDS++c+6Rz7gXn3EsAfgDA/1ABrxi5cqNj68RdC+CT1osJfK01r1GNNxBdBCSYzttHiw/Er1Ffraz/svP0yhildDKsga+9OLjAh0rBkwKeIymUm+PN7k8qySCQAr/G3d8H8DkarndR7J+RaYrROffbAH67WHC1UoHO0U0CfvKeAnZd/B0HPHx/dSGzGaWQYPgulu5+767+cL3GXP2YBg5Yrn9CfazIy0gywh2Av64vd7HlFt2vo/EGkp6MgdUP644+79DipxJ8KYldf0+LPMuuJM2X6Q34SV/CC69yOi5+M0y2C0FZubuuaZfjDWyO663+HBY/lHVyDxhfxxzrL078oTPopbDnRsalgfetCQf4muy8JD6XygV1SOzS8waBlLu+OXZR8IExxLH6FuCPPosMAi3A333eyPVfHHorqz76fEbgY2oFvNS6W0DO+VwyEPh9ig0A1lafAz6QfhJuLL4HYBbjc1Rj/WNaLKbXxCScuGdu4DmLb67rT9+zBl4TX7vjze5PI7++pI1UP6Wxfqy+r1KMH0rzqPGrY1LdoLM5uSPK6ktj/5iWgZ5ka565X7Ql8DFJEneTC69yOm6yMm9m0Lltc9pPzxbE3msDfvjbacCP1+Mn9+aCf3H3PiXpF2oNfE2mPgZ8qJr4nQO7Bm46Stdxl3x74R87FQoM3yF0+VPufitX33fzOYk9TpwvcfX96zgX7+/KC11/oEPoJbDnRsZegA/Vwp1PSQp6DvJSWe4gUEoMbo43UfCBMcCxOD+V4Ns38Hd9e3J9c+AH+Im/bubpJe5KyRXqCXjfyrcGXhNb09Fm91cjbTuxfuZi/el7ZXd/KVd/Wobn6ocarneu+1/iaFHoB9CtYAf2B/j0YhNd/C6x6lJA18cbrJXtc46RGqR6Aj+U5nrZltGBP0ga+8e0iHvvVmTmxvvqFfhQ1u58CXiNFY9BHr53xZyq44YC7njDdvc1cb7U1fcVn5Yru/qS6bxBnGk9aezvqxv3PibJqNYz8D68cwGvteYSqy4tz+nb0hZ/9FnBzd/Wt7P4u89O102tf3+JPKHrkjo5Nx346CIZoUWXwippS+oF+NZ/+G6+1Y9l9yUJvpzF9yVN7AG2Fn/0+ZNrm7ugh5v46wZ6CeylUbBH4GPSTsfVAC8F/eh4fAFdXvB+J+kgQEebidtv6e5nXffKjP62jTbgA3Wuf0yLQ28JO9Av8H55a3c+B7zGmoegpz7jDgCxfsQGgdbg58q3Aj+UFvxdWaH1j2kx6K1hB9oBH6oX4FOw11pzbT3NIBDCbw3+6PPKxJ4GfMlNOoM4A4DU+vtaBnrGKlzpE2NbAi+5gaYW+BrrzoVdC7mm3dJAEIPfMs4PY/ylwAfKN+lcl9cNAFz4F3fvQ0lgT2VATbOsjYC3dOc5sEtBPznO3555fsG7WYrrDayPN02tvu/qhyv3LMAHIF65lwP/us72e0mtf05dQG9h1XefKXatjdVrGcNbZudTwFtDzinPGQhyg0ALq58CHxjDXAv+9j35kl0O+Nt6Muuf02LQW4K+K7OHwFu68xLQpZBr2pQMApbw14APXJe3Ah9ANs73r785BoDFVuRxJHmG3E0FnmvdS8BLIb97lG/v0WX50pF4A0fHjyeuv9blrwE/LK/ZhEOb4BvqDmo1AHTh3oeqhX3bRr/Az2HdNZa8BHquLGcQAPLeQK3VXwp8wC6zP25TNwCU4O8GegvQr9u6ecBzrLu1Na9ti+sNhNZfa/VDdz+M85cAH8gn+Dhz+35bYXvRsgWWFoVe+gho1bPkOgbe0rpzgJdCfvcoPun96JK3EyfXGxj6rrX6EndfCr6vGPiAzb76AFjwD+0BPOsf02Lz9FzgOaADvCm5WHs9A19r3a0g55TjDASlQUBr9aXufmwuX5vVD8sDdZn93Xc0tv6+unHvQ9XAHnu/NAcPLAe8NlkXAt8Kck1bkkHAh19i9TnuPjfObwE+IN+Fx/9s910FAwAH/q6g54IOpGGPfbbvwHPceS7wEtDvHZ9H3394kbkfNXGc3CAggV/q7rcEH+Dfj7+rJ3ygxlBmEGdtf0mLQ28FeqpMr8BbWXcgD7wF5KVyVoPA3aPHRZefY/XnAh9oN5efkjT+j2m5xTnGsMfK3XTgY7BLXXYu6JI2pIOAPwCkrH4p1i/F+TXg+2oF/rYdgUUXuv++llmcw9j3ngt6quy+AM9daFMCngu7BPKnj6Zlv3FZBlo6CNw9uphY/tDqW7j7JfB9hWv1l3iE1u57MQcALvyLu/ehtFZ99BnjqTP7AjzHnS8BXwt66fPSQMAZBIbvILX63Ow+B/zSTTrW++pv3+PcdMODmstOF9DXWvXR5waPi14CeG2yLga8JeSaNiSDQDgAcK1+TZy/j+Bv266P54HFF+fUW/VRuVsEfA3sEtCfuXM2ev3W41Nx+7lBYOizD3/K6mvc/RbgA7b35G/fk6y408fzwFLQr3SZ+GQ5g/gd6AN4jXUvwV4Deelz6SCQGgDuHZ8XrT43u98afCAe59eCf/2ZfACQwF8NPRG9COC/AngXgA2A+865/6BtT+Tqp262qUzYhXVaAN/Suktd9hLokrqlQSDnBaSsfm2cbwV+qJbgX5fhLbqRWH8LS/8YwI875z5PRM8A+BwRve6c+z+SRmqs+vjz8euahN328/bAa5J1NbBLIL/nlX3IsOraQSCEPwf+9j1ZnB9bwTcothvPoHAqD5Btr625NTcmyYq7EkvV0Dvnvgrgq0/+/xYRvQHg3QBY0FvBvi0zfr2PwGvc+RLsUkt+L1E+fN9yEHj66LwIPiCP83354OcW8FjelguUwQfk8G/LdnDDDRG9BOA7APxOrpwjW9i35cavWwIfqkfgtdZcotpBIBwAYuADZXc/F+fnXP3ewAf48G/L6u62M3usFRE9DeBXAHzMOfcw8vmrRPSAiB48Pnu72N7VMbGte81z5TjK7XiTe47c3MA/c+ds98fRvTtnauBz7XHbjPU1NoiF3zs8L3ePHk/On39+w3O/zgzSvnKPAY+WFz5CK/dI89h1nS7LY2WQiaUnoiNsgf+Uc+5XY2Wcc/cB3AeAp55/MTk0aS27RFIrP/4sfSH4F1Buh9pWwHMkgfwdR/HB+W8unxIdp+QBPHPnbGT1Q4sP6OJ8X7FFPIOsMvocSR6aCbSx/NWWnogIwH8C8IZz7qe17dRYdl+9xPG+fEuzBPASC/yOo7d3fzVlpMe3sPjb99Ln2v8dwkHZ//1y3tvU6wuul4K1j7bJ2E1cavlzsnDv3w/ghwF8FxF94cnfh7mVrWAHlovjQ5Wm5ga1BF7qaksgjtWzGACswPeV22Qk5435v2842NeCH10jwnuMgAj+lCyy9/8LrGfWjGXtxtc+THLSniCO17j1rYC3cN+18tuThAC++2/h6uem8ySJPV+tE3tDHaB8Tz4gc/tDLXSXXbmMFnauWsTxoVJufQvgJRadq2fXj6Lvf/3qrug4pQHg3p2zIviAbC6/RXxfkgX4Qz2gHfxm2XtL1QLf2/RcadPK1sBLXPBn1492fzVlpMcPv0PMiwnPAyerP8givi+5+dbaHLVx+7uCntvx3MmwBt5ies6XdKWdFnhprC2BOFbPYgDoCfyU5ojvJ8cUwl9SF9BbwA60BT6URRzfAngN6BrYLdqbC3xf3MReLpu/BPhDO9xMf06LQm8FO9B2AQ5QPz03F/AlScB8bv3N6B9H3OPEBilr8CULd1KyBj/U1Yk9/CktA/3KDnZAfpvs9nN+4k7j1uf3o58XeC3oNWWkx28Bvq/c5qEpNz83uMdUAj9XdtAc8Hfh3oeSwG4BvGQ+nuvW+/IvuDmBl7jaEohj9aQDQEwa8EP551MS3/vSuvklSRbuSOGXqDvoa2AH5ovjQ3Hc+jmAt7bqEknaS/VRCn5pAQ83vl/KzS/NAHDhl1j9bqDnxu25E1A6gVLga6fnci5lK+BLkoD53Prt6B9H3OO0At9XLr73xXXza8HXyBL+xaGvids5bSc/EwAfijs958u/8OYGXgt6TZnY8VNqAT43scd90q8U/Jxi1p475y9x+1NaDPrauH3UVuSkWSbuaqfnlgLeEnSLurnBJ+buW4OfUsrN5+ZsYtJk86Xwa7UI9KxluFx3prC/3baMTeIulCaO99UCeI1Vt5KF9W8J/r64+Vz4tVZ/cfc+lAR2C+AliTvtMttB/gXYCviSJJb576zejv5ZHmMO8H3N7ebn4C15uq3g7wb6WtjjZfWJu1CaZbYpt35u4LWg15SRHJszWNWAb+3ml8DPSbpoZ1uHDz9Hi0NvCbtkvlLzCOmYNHG8r5bAW4JeqstRbgAIvwNnJkIyj69x831J4nvr1XqpdmLi8LTsMlxDy94yjq+dnvPlW6AWwGusuoWkA0esjyXwpWv1a938nEdn6eYD/MS2xNNNaRnoSZ+ki5YzjONDaTK4XLe+JA3wHIngXJ2P/szbNwA/J42b70vi5kuUSupZWv2UFnfvU9ICXyvtXndS1TxVpkYSyx6DnAu+9FiWqj232t+2dpnurl5j8LuFniMN8JLkHVdc1z6n0E1tYeVrged8VnPMXftCa2+xjbfFb1hS6003uOoSesuTU9OWpWu/tKyAl5ThHttyrYBEpU01Y6pZsBMqu2CsobXvDvpat75mnTPXtddIksCz3rDSGnhN2daufs7a18b1OdVM35XUCvzuoOfIKo7X/kC18bxUNa59K+A1dUR9qUzoWcb1NQO+Nq5vqa6gt3owQLp9+x+Au+R2CbUGXlM31aelXHxfPcb1Lax9N9BbxPGSufmSamM3rhspTeCFSsEyF/CaNqxcfW0YxL3fPqU5pu5GZYxnqLqAfqnpOV9LxfNScZaszg28pq1YH6Vz9qG4cX1rWcb1EnE56gJ6jmqBt5qqmzuel2op4DVttk7uWa6FWDqutzR4i0Pfau5yDte+Np63nptvAfyzq+s/67Y1stquOyXNnXdStdhZJ9V2TItCb+nWt5qq06hm6e0ckgCfe21xjHCgKrn4obQLdXqK67mySuotbulLahHHS1z72nje11JLb31pgS+9X3OsudRzXM/eYMOAh+W2yzKcnpvrBhvNXXUpLeHa1wLP/VxzzFbqYaAF+pqv73a7LK5q4yELt6zXpbe+rICXluMcW+riz7EWf1BukY7lHnrb9+ax9t26963i+DmnUyyn6nxJrbw18JryS1v8kjQba4SyMCCWSb2UuoTeMo6vcass43lfta69RK2A19TrAfyW6/BDtTQw+/csOwMlNyFoNFVnGc/PpdbAa+qn+tTSxd+3+frW1r476FtPz/lqGc9zM8WtEnhzAa9ppweLbyXLW22l0lp7k5+ciD5ERH9IRF8mop/QttNymS3Qx1SdZZIppbmB17TXO/jaRToSQ5LdHrthUq/6ZyeiNYCfB/C9AN4L4AeJ6L3SdrTTc+PPwtf2U3W9qwXwz674m6rXgG/t4nNUu0inpJq4vpWbbzHWvw/Al51zf+ycuwDwSwC+36DdiSxPgta1t7yV1tq1bwl8K/DnEneXXI5aJXi1WuJR1e8G8Gfe6zefvDcSEb1KRA+I6MHVN5e/d/qgm685Qql9lAX0saU2E5PsnLvvnHvZOffy+lueUh1odWG3qsdd1n/1R5d1dwv9zeX4PHz96u7o9d9efUvwOn/e/nbDs8ZfFzg5X9+cj/61br+lHj4+3f3/Le//37ise+zr5cUddlm6aO/2rITRp0WP3gTwovf6BQB/Lm1E2vFoG402rjm/uPafHl3yf/DW+tpmOgi0BN+63bCv4fcJB7lwEAwHyXAQ5ejhxXUfagdxwMaYDOIYOQ03Fj38XQD/iIjeQ0THAH4AwH8zaDcqibVfCUbZqxlGZGBsfZaWtUWuAb43+YO7P+iHqrluckbK0qudtF3bgHPuMYAfA/BpAG8A+GXn3JdUnWmcJJe4WhIXrqS3MqBLXXyOJEBZgd8b8K0G15rrgmuEuMBreTExb865/+6c+8fOuX/gnPv3Fm3m1HIUrJHvKtbGjSmFLm/MxQfmBd8C+JaufSqe93+vFmoZz9cYyO4mV7hfJgW+JK7Xxl++62cRB7bSHOD3ZuG18n9Hbd5Gcj2lrtNWcfyofl31NrJ082viet+Vy8V1UoWupyYB5Stl7YG24C8FfGjlW8oqnudch3N5sF1Cz1Wvbn5Mubg+VO3UXagW4FsCL3XtQ4WDZmqqrlY9TNWZzHLVN9FGrZJ6c8ybDmoV12tkCf6+uvSpeL6HqTpLt74U4i4CPRkuKY6drJq4XjMF41803ORQbXaZm9Ab1TEA3xp46wRe7XnlxvOWU3WWHivn2u/W0gP1ST1rtVyk02LqLqYa8PfVwmvl/941rr0kr5QShwWusVsM+lar5+LH0n/Nuebr55QU/OHPun2Od1Ijy6W3OVmuwovJOtTt2tIDbax9OBJb/2it5utLCT0JRK0scQ3wc7n2tfH8nK69VRw/Kssvai+2OzLz7ey1cT1Xpak7zdTUkuD3YuG1anFfhYVrXz6GsHybbgg60MjNnyN86OnmG19LgN8b8LWhlFU8X6NWxm5x6LmqXoUU/Bjadfi1i3SkF6N2zn5O8C2Ab+nat1x6WxMa5lz7Fm79ro68ir16cPOXiuutV+f5mgP83iy8REvH87XSttcF9FaynLqbK663UG6FXkvwrYCXWvkW4t5Km5PVVF3rHFY30Lew9tNMqf3XbRHXW8/ZtwB/SQtfcu3n2iWnxVRdS7d+V1df1V5zzt0D88b1/sUXxvVzbKxhCf6+uvS18Xzrqbo5gAc6g56rVu6PxZJcXy3v19bchCMFPwa3NfBzbIkVk8WttDm1mqqzMIzdQV/r5ocj6NxTd9y4vmTtLebsY/ra5imV1U8NAqljzKWca68Rd6ouNBC1U3WWy2zXhS0Nu4Me2H451o0Dl4qdQAtTdzlrr3HxJavzpODHrD33tlsJ/NwbZ7htDv20nqYLgefcVWdxg00pgVdy7a2AX5+XgQeWgt4xO6e0+iVrXzNn7yt1A07O2udie6Ae/O17cahiksCaq19SqU9LAB8qZeVD4HMJPGvgOQaQC/uuTX5Re3HB51r9cT0Z+L641j4n7u22HHe0JpvfYgDQWPV8mXbA55SappMAnzMYNcBzr3sJ7Lu25VVsxXZJGCehBH5OkptwUm5+zlXMufmabL5mGk/j/g9gW7jv8bJtgee49dxQrQR8zpBwgZfArgEe6AD6QZbw5+qOX/Pj+1xMV+vmA/LEHhAHXwK/ZACwbHPoZw/Ah0pZeSnw2T3tI8DPAfvuWHXV7cX9QqmTZB3fp8C3cPNbgA+koYpJOgCk6pfL5fvUEvicuG69r9rEXSjTJB0jud0d9IBsNIvBX+PmS1Tr5gPtwB/EhX9bljcAaKx6Tq1jeI5br115V5u4s0rSSWayFtsjj9NBKfyj14Ukyfi1jZvvK+fmzw0+ILP+2/JjsC3c95haAF8bx0vcel+WwLeAfVdHVtxW3A6bxDGG03jcufslwOcs4tEMAJZtDv2cE/icrNz68Wd1wJekgX1XV1fNVtzOl+CXuvkS8DXZ/FAtwK+BH5APAKn6JaVAH2QNfKja6bmaxB1XHONWA/uujbrqdpJ8mdzJkbj50fJM8DXZ/O3rvPXRTOel1qOXQAvFHQA0Vj2nFsDvUxwvceU5WhcGoO72exq+2IYxdTqcqCvBfS2rC8Lm2Gbj/auLFdbH221iLy/u4Oj4cbTco8s7uHsU/+zhxQnuHed/8bcen+KZO2fXdR6f4p73GtiC846jt5Nt+OA9u36UPR4wnj9/bv1NkScguU9gTuBz4sbxLYAvyQr2XXu8YsZibKcstfyjuh26+TXxPVBn8UNpPADLNoewpDXwoSym53y1Bp6d87rgAw8s6N6zRyVBsm9Uzxh8Xxo3P1RL8FOxfkzSAaCmbq5fLYCvdetrEnc1agX7rn15FTtJOs05ES3B12TzAVl8HwOfE+Onlu22GgA0oOf6MCfwoXJuvS9p4k5j5VvDPqiLmH74AleMGRZJzM9RGOOvLoDNsf/5CpsncTtdrOCe/N9drkBH13GKNr5/dHmMu0fXv2Asxv/G5QmePrp+L4zxgXic78uHLhf7D6q5f1+y0UVr4ENp3PrWcbx1zF5SN9l7wMbyS619jSzi++3rvMUHyll9gL/tltQDsG5z8E7mAN7Crfe1BPASLtYXDuuLfKK6Cnoi+iki+gMi+n0i+jUieq6mvUG18M/p5ufi+6XAz7n8oWoGAGndXL/mBl7j1lsvseXIEvZBtZb+dQDf6pz7dgB/BOCTnErcR1VL4R/V7SC+DzUX+Lv6jQYADei5PrQGPpR2eq5WEivPvfYlsA+qiumdc7/lvfwsgH/JrTt09Oq4nPWUxPw5TeL1y3FuoBTf+9LG9+cXRzg5vv61OTE+gFGcz4nxQ/nQ5WL/QTWuv2R339YxPNDn9FwKeIll18pyKPtRAL8prSQZqUqjX8nac+rkLL7F/H1MnFV7nCW73E0hpR5AizZj/bUC3np6bg7gW1r2UEVLT0SfAfCuyEevOed+/UmZ1wA8BvCpTDuvAngVAI6fesfkcyvLH1rv9fl4xV7Oel+XSa/a87P5oXyL71t7YGzxQ2sPlC0+MM3shxYfGFvOkvUHplaZ4wWk6uZUGpA498NLXfoe4vgS8NaWfX1eLleE3jn3wdznRPQjAF4B8AHnXPKIzrn7AO4DwFPPv5gsZwG/FPyw/LbMNfjcaTxABj6AoqsPIDulN8ASwg9MQZMOArEBwBL0QSXgY56PJfCpOL4H4C1hH1QV0xPRhwD8OwDf6Zz7Jrseo3/+ly0NADH49wF8oBzjb9+Tz+XHJB0EpO6/dM/5Fu48sP/AtwDdV+3inJ8DcALgdSICgM865/4Np6Lf4auTEtQ867++uDngA1Al+AaVBgBA5wmU2sgp+/ReA3ce2G/gW8M+qDZ7/3la83IAAA6fSURBVA+rjv5Ew5fgwg+kBwAL8IF0Vp+b0Q/FAR+od/cHhYBZDQJWkPtq4c4D+wM8B3YJ6KWyXSzDHaSx/sB0AKgFP1Zn9JlXPkzscafygPhyXa27DyB7i67FIMCR5CmxGusO3C7gubDPFtO3FNf6A3H33xr83Bx+LfgAqt19YAxR6R59zSDAaSen0g611u48sL/At4B90DLQbwQuvRD+XsEHkIzxgXKcD/Dc/UGSAQDgDwLS571ztqJu4c4D/QKfgr0l6L4Wt/RS+EtlLcAH0gk+LfjANLkHoDrO3743tfq+QvA0gwBX0sdzt3Dngf0Cfi7YBy0O/SBRPF8oWws+kM/sW4EP2Lr7g3IDAKAbBLhtlZTb2cbCnQduFvCiBB4z+98N9L5ELn2ibAvwx5/xV+1JwQf4Vh9A1PL7kg4CQH4gkIDO3Ybawp0Hbg7wLWAftAj05NwEypik1t8C/O2xvDK+Fc/coJObwwd44AN5dx/gxfqhJF7AIKkFjx2LV76NdQfaAJ/a8cZXDfBsV79i/f2ill605JYxAFiAv21HDz6gd/UBvrsP8Kx+qBiU3IFA0ma+fPyys7LuwDzAp26cGb8uAz8H6L66cO8lS26BvPsfA396vHbgT8o+ueByK/cAiN19IG31fZUGgW0dWThgBbkvi2Qd0O5uOYspuRjcHOAlsK/Py1tNdwG9LwvrH4LPnccH+FN60nvxwyW7AFTuPlC2+r4sBgGpuHvNA+3ceWC/gZdadQ7sg5aB3viGGyACusLVB2Rz+TXgAzp3H+Bb/Zg0g4C0TY44sAPyZN2gfQXe2qrHtJilHzp8dVLex4M7AMwFfq6sBfjA1N0H+FZ/kGYQALjegOzSyW1fBdhZd2A/geffbKMD3dfi7r3/JSwGgBj427bH7n4I/ra963Zy4E8SdgbgA1N3H+DBn9puWzMIADrLXTp2Ti3c+e1rBK/7A74F7Ouzq+zni0Pvy2oAiFr4SJxveaNOmBOQgg/E79LTZvhjioHIHQhKkkAOpLex0iTrgGWB507JSYGXWvUS7IO6gt6XZgBYGvywPAd8AKZWf1BpABik9QasIB9k5c5vX4efLwu81J23tOoxLbM4Z+OwPrvC1emaVZ4b/8fA39arz+zXgA8gez++pdUfFIOMMxCkvAEp5Kk+xGSRrLt+L3y9H8C3Bt3XsotzvM5zBgCO9Y9ab4MEX+0uPLlFPNv3bK1+TCGEWm+A235JltZ9+174un/g54R9UDfuveUAsA/gx8pv36uz+kAZ/kHaQSBVv6TcnvOAnTu/fe9mAC8FfXVe/g27gd7X8EVr3P/bCj4whctyELAGfZCVO799b/+BbwH7oG4X5wBy6z+p3wn4oazBBxCFf5DVIMAVF3TA1p3fvqd/vlwPwEtgl4A+qqeqZaDV+WNRp9dnV8UTEjuZZj9UGC8KnqYTuzEj9iSd2PPzYk/ViUGRe3hDqKuL1e7PStI23eXKPH6vuT12aeA51/cgKTuhFnfvh85vTnhdKbn+6/MNK8af1rO5Q2/3WcHN35aZ3pNvYfWBvOX3pfUCpAMGZ1AqLbbZvjetF8IOyIAPNSfwc1j2UItDP8j/QpwBIOf6c8DnzONH60WW7V7Xl8X32zK24A+axMgGg0AL0AFbdx6QA++X7xF4Keyrs3z5xdz7nCxcf46rb7EuuvYR2dsykQs3eoHH3eHc01YHDe60JAwA9G77TQA+prmAHxjgcrA6e7z7K2mxxTkcaay/b/UtLf64jfqMPjDN6gOYJPi27/nlpvP5QN7qh9J6Adz2Sso/Elo+HQekBtPxa8nz4YGykWgBvLVVj2kx937o7OaU1wVp7D/ICvzW22tfl6tz9321GgSkkAN50AFb6x6rWzs11xp4Cewa0H0tHtNr4QfiA0Bsee9NBh+YWv1BPmjcAQDQQZ07dkm3Gfg5Yd+1Y9KKgSQxya5OIuaJnVxujF/6gWun8qIXsyDO58b6oYbYn5sD0EjT/gH4sqRcFNsza8lQFvBrwd+W6xf8ZP3InH5OVoOAFvRYf+NrFeLz77cBeCkHdH65+8tpoRV5zESe94U57v/q/PHI5de6+tty+ZV7rZ6kA5QTfEP92Ko/H6SU2x8TNx+gGSA4g5GldQf2F3iNRS9BPjmG+AhG4oxIvrij3j5ZfKCNu3/9+SppVUsKPQErax4vG3vvdgHfyqpHjyWuYSxp5zWx/00CnzunH1PNINCibWt3HrAHPqYWwHOlBX10vKraxrKy/pbJvW25wrZGDcG3svox1Q4CmrpDP2ut+1zAFwf5mYC3gH3Q4lN2MQ1fzp3wNm+ITfuF8T1QF+OP67R7hBbAn9bblp3G+r5yd/uFisHr5wV0AwO3nJ07H6vLfVz0uN1lgVdD/ugs+/Eyln7DSzBpXP/Ra+PpvNzraB3GEtwWVj9sK2ddS7K25vHyywMvDuMaAq+y6o/Orv8KMoGeiD5ORI6InmdXYnZwdwzmibhJ4NfG+jHVDgKWbVrH70P9UZk9AV6dmBNyBBi490T0IoDvBvCnqgb8Dt89LR/POykp93919nji6gMwmc5rsQnHtp34AzVSq/gAvsufUwwmTkig8xoKXgkTdkAHfLydPoAXSQj5pA9Vtbf6GQCfAHs/nIwMrT8nwZe6Oy/8Yee0+Fx3P9XGtjyN/qSKWe5aa24FfPIcGTyMonvgFVY92o+aykT0EQBfcc59kVH2VSJ6QEQPLq4e5QsL4hMgDb9lZr8V+C3c/Wnd+kGAX1Z2LAt3nnMO9hZ4IQscFd17IvoMgHdFPnoNwE8C+B7OgZxz9wHcB4BnT/4e3yt4dMZy+5cQx9XXt51+fl6o3A4+0fKKMEDSbn070/dyO91Myioy9S1VBbxGZ/mTVYTeOffB2PtE9G0A3gPgi0QEAC8A+DwRvc859xfZRpn30+/EBJ/OLydxfhjfA/HpPI44226V2yg/MDOn0qabElnBr4WdY6Gl9VvE8dE2jPagT6oR8ECFe++c+9/Oub/rnHvJOfcSgDcB/NMi8ILOjWTo3sTEzegX21E+oDDeFr9sDUA1FloVLmRc8lorP6lrYOW510HMtVdZ+YbAA0vN0w/SgN8Yfqk4yzY5kkArTezx2pTH+1rgLdTKyi+uxsADhtA/sfh/La54dm5q9edO6GnEuTivjy9rm5vYStbnJN80AwSjX9ZWXiOta29i5WcAHlja0vuSgt+xtC6+hbXXtDVtOw11K+sumiFgWnnNzTSTNhQhnkozAQ/0BD0g+wIzWvuSVBeTgbXngG8Fv3q6r+L4nPn4Glkl8LhWPqkZgQd6gx4wA99SFi5+C2vPVXWW3DhZN21f3PxIHCuvatfYylvdJQegyjPuD3rABHzTE6wQ775sSXvx99k3tVRafYksjmNt5S0SeLUbW0alMVyVoXCf0AOyBB/zxFm6+BabL6SUurgtklot4de03cLK69qxy+KbrK9PySD31S/0gzpJ8FnM2UfLGEy/qW5+MQZf017NzTQ57VUCTyIjFrq+n16kiLXfh4Se/Bjx97XgW8DfGviUtFZ+LxN4hsZvMUvvFvrClrJK6FlY+xpp4VfXk27rVWHlNeougWd8/S/q3puDz7T2WrVy8YF5k3opiTLuSthzfWxh5edM4DWx8g0M3uIxvXt0xoffKqZZ2MWXWCNpUs9iR5wc0HNZd05fQrWcStWqyug08nAXh37Q3OBbSbPBYrqt2t5cq4XVb2HdB9Vaec0uxvF+lL256mk6jhpe591ADxiB3zihx7kolrb2u3pGVn9u6+4fO1TNwGjl2sdk6to3NmxdQb+vsnQRLa39IMtNMCXH5Kp2/YH2/GvqSay8yrWfwZPtDvoerX0oq4TeHNZ+1MYM8Fsdo9bKaxJ42qy9mZWfKXTtDnqgv/h+7jn71lN4rcDXtGu9EMeyXs3uOGIrb3gtl/jpEnqgP/BL0ib0bOaV5XVa7HsvVauFON0m8FLX9IzAAwA55mOjLUVEfwXgTxo0/TwA+UYey2if+grsV3/3qa9Am/7+fefcO2MfLAJ9KxHRA+fcy0v3g6N96iuwX/3dp74C8/e3W/f+oIMOaqMD9AcddMt006C/v3QHBNqnvgL71d996iswc39vVEx/0EEHlXXTLP1BBx1U0AH6gw66ZbqR0BPRx4nIEdHzS/clJyL6KSL6AyL6fSL6NSJ6buk+hSKiDxHRHxLRl4noJ5buT05E9CIR/U8ieoOIvkREH126TyUR0ZqIfo+IfmOuY9446InoRQDfDeBPl+4LQ68D+Fbn3LcD+CMAn1y4PyMR0RrAzwP4XgDvBfCDRPTeZXuV1WMAP+6c+ycA/hmAf9t5fwHgowDemPOANw56AD8D4BMAus9QOud+yzk3rO38LLZP/u1J7wPwZefcHzvnLgD8EoDvX7hPSTnnvuqc+/yT/7+FLUzvXrZXaRHRCwC+D8AvzHncGwU9EX0EwFecc19cui8K/SiA31y6E4HeDeDPvNdvomOIfBHRSwC+A8DvLNuTrH4WWwM169a78oe0Lywi+gyAd0U+eg3ATwL4nnl7lFeuv865X39S5jVsXdNPzdk3hmKPtunegyKipwH8CoCPOeceLt2fmIjoFQB/6Zz7HBH9izmPvXfQO+c+GHufiL4NwHsAfJGIgK2r/Hkiep9z7i9m7OJIqf4OIqIfAfAKgA+4/hZNvAngRe/1CwD+fKG+sERER9gC/ynn3K8u3Z+M3g/gI0T0YQCnAO4R0S86536o9YFv7OIcIvq/AF5WPT57JhHRhwD8NIDvdM791dL9CUVEd7BNMH4AwFcA/C6Af+2c+9KiHUuItqP9fwHw/5xzH1u6P1w9sfQfd869MsfxblRMv4f6OQDPAHidiL5ARP9x6Q75epJk/DEAn8Y2KfbLvQL/RO8H8MMAvuvJ+fzCE0t6kKcba+kPOuiguA6W/qCDbpkO0B900C3TAfqDDrplOkB/0EG3TAfoDzrolukA/UEH3TIdoD/ooFum/w+AWKL9UhmTmwAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
"plot(T)\n",
"print(\"Temperature at center at t=0: \", T_0(0,0))\n",
- "print(\"Temperature at center at t=%.2f: \"%(k*delta_t), T(0,0))\n"
+ "print(\"Temperature at center at t=%.2f: \"%(k*delta_t), T(0,0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See this [tutorial](https://fenicsproject.org/pub/tutorial/html/._ftut1006.html) for a more extended version of this example."
]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "jupyter": {
+ "outputs_hidden": false
+ },
+ "pycharm": {
+ "name": "#%%\n"
+ }
+ },
+ "source": [
+ "The [`solve`](https://fenicsproject.org/olddocs/dolfin/1.3.0/python/programmers-reference/fem/solving/solve.html) command is a high-level call to solve a linear variation problem `a==L` or a linear system $Ax=b$.\n",
+ "\n",
+ "Alternative ways to call the solver for *linear* problems, including manual assembly of the linear system, are discussed in this [tutorial](https://fenicsproject.org/pub/tutorial/html/._ftut1018.html#___sec121).\n",
+ "Solution of *non-linear* problems is briefly introduced [here](https://fenicsproject.org/pub/tutorial/html/._ftut1007.html#ftut1:gallery:nonlinearpoisson) and discussed in more detail [here](http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/nonlinear.html)."
+ ]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"pycharm": {
"stem_cell": {
"cell_type": "raw",
"metadata": {
"collapsed": false
},
"source": []
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
diff --git a/FEMII/01_fenics_misc_topics.ipynb b/FEMII/01_fenics_misc_topics.ipynb
new file mode 100644
index 0000000..42587a5
--- /dev/null
+++ b/FEMII/01_fenics_misc_topics.ipynb
@@ -0,0 +1,60 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# FEniCS How-Tos"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Subdomains (cells, facets, edges)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Boundary Conditions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## File IO"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/FEMII/02_hyperelasticity.ipynb b/FEMII/02_hyperelasticity.ipynb
new file mode 100644
index 0000000..44fc0bb
--- /dev/null
+++ b/FEMII/02_hyperelasticity.ipynb
@@ -0,0 +1,293 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "pycharm": {
+ "name": "#%% md\n"
+ }
+ },
+ "source": [
+ "# Hyperelastic Material model in FEniCS\n",
+ "\n",
+ "Based on the FEniCS [hyperelasticity example](https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/hyperelasticity/python/documentation.html).\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from fenics import *\n",
+ "\n",
+ "# Optimization options for the form compiler\n",
+ "parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
+ "# ffc_options = {\"optimize\": True, \\\n",
+ "# \"eliminate_zeros\": True, \\\n",
+ "# %\"precompute_basis_const\": True, \\\n",
+ "# \"precompute_ip_const\": True}\n",
+ "\n",
+ "# Create mesh and define function space\n",
+ "mesh = UnitCubeMesh(24, 16, 16)\n",
+ "V = VectorFunctionSpace(mesh, \"Lagrange\", 1)\n",
+ "W = TensorFunctionSpace(mesh, 'DG', 0)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Mark boundary subdomians\n",
+ "\n",
+ "left = CompiledSubDomain(\"near(x[0], side) && on_boundary\", side = 0.0)\n",
+ "right = CompiledSubDomain(\"near(x[0], side) && on_boundary\", side = 1.0)\n",
+ "\n",
+ "# Define Dirichlet boundary (x = 0 or x = 1)\n",
+ "c = Expression((\"0.0\", \"0.0\", \"0.0\"), element=V.ufl_element())\n",
+ "r = Expression((\"scale*0.0\",\n",
+ " \"scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])\",\n",
+ " \"scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])\"),\n",
+ " scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3,\n",
+ " element=V.ufl_element())\n",
+ "\n",
+ "bcl = DirichletBC(V, c, left)\n",
+ "bcr = DirichletBC(V, r, right)\n",
+ "bcs = [bcl, bcr]\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define functions\n",
+ "du = TrialFunction(V) # Incremental displacement\n",
+ "v = TestFunction(V) # Test function\n",
+ "u = Function(V) # Displacement from previous iteration\n",
+ "B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume\n",
+ "T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Kinematics\n",
+ "d = u.geometric_dimension()\n",
+ "I = Identity(d) # Identity tensor\n",
+ "F = I + grad(u) # Deformation gradient\n",
+ "C = F.T*F # Right Cauchy-Green tensor\n",
+ "\n",
+ "# Invariants of deformation tensors\n",
+ "Ic = tr(C)\n",
+ "J = det(F)\n",
+ "\n",
+ "# Elasticity parameters\n",
+ "E, nu = 10.0, 0.3\n",
+ "mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))\n",
+ "\n",
+ "# Stored strain energy density (compressible neo-Hookean model)\n",
+ "psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2\n",
+ "\n",
+ "# Total potential energy\n",
+ "Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds\n",
+ "\n",
+ "# Compute first variation of Pi (directional derivative about u in the direction of v)\n",
+ "Fpi = derivative(Pi, u, v)\n",
+ "\n",
+ "# Compute Jacobian of F\n",
+ "Jac = derivative(Fpi, u, du)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Solve variational problem\n",
+ "solve(Fpi == 0, u, bcs, J=Jac)\n",
+ "\n",
+ "# Save solution in VTK format\n",
+ "file = File(\"displacement.pvd\");\n",
+ "file << u;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "F = variable(F)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "S=diff(psi, F)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "VariableDerivative(Sum(Product(Sum(IntValue(-3), Product(Power(Determinant(Sum(Identity(3), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 107), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 118)))), FloatValue(-0.6666666666666666)), Trace(Variable(ComponentTensor(IndexSum(Product(Indexed(Sum(Identity(3), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 107), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 118))), MultiIndex((Index(99), Index(98)))), Indexed(Transposed(Sum(Identity(3), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 107), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 118)))), MultiIndex((Index(97), Index(99))))), MultiIndex((Index(99),))), MultiIndex((Index(97), Index(98)))), Label(3))))), Division(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 123), IntValue(2))), Product(Division(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 124), IntValue(2)), Power(Sum(IntValue(-1), Determinant(Sum(Identity(3), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 107), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 118))))), IntValue(2)))), Variable(Sum(Identity(3), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 107), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 118))), Label(4)))"
+ ]
+ },
+ "execution_count": 31,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "S"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Expecting a Variable or SpatialCoordinate in diff.\n"
+ ]
+ },
+ {
+ "ename": "UFLException",
+ "evalue": "Expecting a Variable or SpatialCoordinate in diff.",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mUFLException\u001b[0m Traceback (most recent call last)",
+ "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0mS\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdiff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpsi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0mS_fun\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdiff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpsi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mF\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/ufl/operators.py\u001b[0m in \u001b[0;36mdiff\u001b[0;34m(f, v)\u001b[0m\n\u001b[1;32m 359\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mVariableDerivative\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 360\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 361\u001b[0;31m \u001b[0merror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Expecting a Variable or SpatialCoordinate in diff.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 362\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 363\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/ufl/log.py\u001b[0m in \u001b[0;36merror\u001b[0;34m(self, *message)\u001b[0m\n\u001b[1;32m 170\u001b[0m \u001b[0;34m\"Write error message and raise an exception.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 171\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_log\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 172\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_format_raw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 173\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 174\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbegin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;31mUFLException\u001b[0m: Expecting a Variable or SpatialCoordinate in diff."
+ ]
+ }
+ ],
+ "source": [
+ "def Energyfunction(C):\n",
+ " I = Identity(3) # Identity tensor dimension 3\n",
+ " F = I + grad(u) # Deformation gradient\n",
+ " C = variable(F.T*F) # Right Cauchy-Green tensor\n",
+ " J = det(F)\n",
+ " I1 = tr(C)\n",
+ " I2 = 0.5*(tr(C)**2-tr(C*C))\n",
+ " I1bar = J**(-2.0/3.0)*I1\n",
+ " psiV = (lmbda/2)*(J - 1)**2\n",
+ " psiD = (mu/2)*(I1bar - 3)\n",
+ " psi = psiV+psiD\n",
+ " return (psi, C)\n",
+ "\n",
+ "psi, C=Energyfunction(u)\n",
+ "\n",
+ "S = 2*diff(psi,C)\n",
+ "S_fun = diff(psi, F)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Stress: 5.511622533084147\n"
+ ]
+ }
+ ],
+ "source": [
+ "print('Stress: ', S_fun.vector().get_local().max())\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Coefficient and argument shapes do not match!\n"
+ ]
+ },
+ {
+ "ename": "UFLException",
+ "evalue": "Coefficient and argument shapes do not match!",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mUFLException\u001b[0m Traceback (most recent call last)",
+ "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mS\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mderivative\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mPi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mF\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/dolfin/fem/formmanipulations.py\u001b[0m in \u001b[0;36mderivative\u001b[0;34m(form, u, du, coefficient_derivatives)\u001b[0m\n\u001b[1;32m 80\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Computing derivative of form w.r.t. '{}'. Supply Function as a Coefficient\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 81\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mderivative\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mform\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcoefficient_derivatives\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 83\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 84\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/ufl/formoperators.py\u001b[0m in \u001b[0;36mderivative\u001b[0;34m(form, coefficient, argument, coefficient_derivatives)\u001b[0m\n\u001b[1;32m 279\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 280\u001b[0m coefficients, arguments = _handle_derivative_arguments(form, coefficient,\n\u001b[0;32m--> 281\u001b[0;31m argument)\n\u001b[0m\u001b[1;32m 282\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 283\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcoefficient_derivatives\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/ufl/formoperators.py\u001b[0m in \u001b[0;36m_handle_derivative_arguments\u001b[0;34m(form, coefficient, argument)\u001b[0m\n\u001b[1;32m 223\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoefficients\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marguments\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 224\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufl_shape\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufl_shape\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 225\u001b[0;31m \u001b[0merror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Coefficient and argument shapes do not match!\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 226\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mCoefficient\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mSpatialCoordinate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 227\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;32m~/anaconda3/envs/fenicsenv/lib/python3.7/site-packages/ufl/log.py\u001b[0m in \u001b[0;36merror\u001b[0;34m(self, *message)\u001b[0m\n\u001b[1;32m 170\u001b[0m \u001b[0;34m\"Write error message and raise an exception.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 171\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_log\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merror\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 172\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_exception_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_format_raw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 173\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 174\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbegin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+ "\u001b[0;31mUFLException\u001b[0m: Coefficient and argument shapes do not match!"
+ ]
+ }
+ ],
+ "source": [
+ "S = derivative(Pi, F, u)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.3"
+ },
+ "pycharm": {
+ "stem_cell": {
+ "cell_type": "raw",
+ "metadata": {
+ "collapsed": false
+ },
+ "source": []
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/FEMII/hyperelastic_2.ipynb b/FEMII/hyperelastic_2.ipynb
new file mode 100644
index 0000000..afc3f87
--- /dev/null
+++ b/FEMII/hyperelastic_2.ipynb
@@ -0,0 +1,493 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Hyperelasticity\n",
+ "\n",
+ "inspired by https://gist.github.com/alen12345/da8e4bcabe6eb7c6659ee9ef8a323a32\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from __future__ import print_function\n",
+ "import subprocess\n",
+ "import fenics as fe\n",
+ "\n",
+ "fe.parameters['linear_algebra_backend'] = 'PETSc'\n",
+ "fe.parameters['form_compiler']['cpp_optimize_flags'] = '-O3'\n",
+ "fe.parameters['form_compiler']['quadrature_degree'] = 3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ " \n",
+ " FEniCS/DOLFIN X3DOM plot \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "mesh = fe.BoxMesh(fe.Point(0,0,0),fe.Point(10,1,1),100,10,10)\n",
+ "mesh"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "V = fe.VectorFunctionSpace(mesh, \"Lagrange\", 1)\n",
+ "W = fe.TensorFunctionSpace(mesh, 'DG', 0)\n",
+ "\n",
+ "# Finite element functions\n",
+ "du = fe.TrialFunction(V)\n",
+ "v = fe.TestFunction(V)\n",
+ "u = fe.Function(V)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Mark boundary subdomians\n",
+ "\n",
+ "left = fe.CompiledSubDomain(\"near(x[0], side) && on_boundary\", side = 0.0)\n",
+ "\n",
+ "# Define Dirichlet boundary (x = 0 or x = 1)\n",
+ "c = fe.Expression((\"0.0\", \"0.0\", \"0.0\"), element=V.ufl_element())\n",
+ "bcl = fe.DirichletBC(V, c, left)\n",
+ "bcs = [bcl]\n",
+ "\n",
+ "# surface load\n",
+ "class RightBorder(fe.SubDomain):\n",
+ " def inside(self, x, on_boundary):\n",
+ " return x[0] > 9.9-fe.DOLFIN_EPS and on_boundary\n",
+ "\n",
+ "boundaries = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim() - 1)\n",
+ "boundaries.set_all(0)\n",
+ "right_border = RightBorder()\n",
+ "right_border.mark(boundaries, 1)\n",
+ "boundary_file = fe.File(\"boundaries.pvd\")\n",
+ "boundary_file << boundaries\n",
+ "\n",
+ "dsp = fe.Measure('ds', domain=mesh, subdomain_data=boundaries)\n",
+ "load_s = fe.Constant((100., 0., 0.))\n",
+ "integral_S = fe.inner(load_s, u) * dsp(1)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "########################## VARIATIONAL FORMULATION ###########################\n",
+ "# Large displacements kinematics\n",
+ "def largeKinematics(u):\n",
+ " d = u.geometric_dimension()\n",
+ " I = fe.Identity(d)\n",
+ " Fgrad = I + fe.grad(u)\n",
+ " C = Fgrad.T*Fgrad\n",
+ " E = 0.5*(C - I)\n",
+ " return E\n",
+ "\n",
+ "E = largeKinematics(u)\n",
+ "E = fe.variable(E)\n",
+ "\n",
+ "# Stored strain energy density for a generic model\n",
+ "def strainDensityFunction(E, Ey, nu):\n",
+ " mu = Ey / (2.*(1+nu))\n",
+ " lambda_ = Ey*nu / ((1+nu)*(1-2*nu))\n",
+ " return lambda_/2.*(fe.tr(E))**2. + mu*fe.tr(E*E)\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Using def gradient, neo-hookean\n",
+ "\n",
+ "def defGrad(u):\n",
+ " d = u.geometric_dimension()\n",
+ " I = fe.Identity(d)\n",
+ " F = I + fe.grad(u) # Deformation gradient\n",
+ " return F\n",
+ "\n",
+ "F_el = defGrad(u)\n",
+ "\n",
+ "def strainDensityFunctionNeoHookean(F, Ey, nu):\n",
+ " mu = Ey / (2.*(1+nu))\n",
+ " lmbda = Ey*nu / ((1+nu)*(1-2*nu))\n",
+ " C = F.T*F \n",
+ " I1 = fe.tr(C)\n",
+ " J = fe.det(F)\n",
+ " return (mu/2)*(I1 - 3) - mu*fe.ln(J) + (lmbda/2)*(fe.ln(J))**2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# growth\n",
+ "\n",
+ "def defGradGrowth(stretch_ratio=1):\n",
+ " return fe.Identity(3) * stretch_ratio\n",
+ "\n",
+ "F_gr = defGradGrowth(1) \n",
+ "\n",
+ "F = F_el * F_gr\n",
+ "\n",
+ "F = fe.variable(F) # this is needed for stress computation, to be able to differentiate energydensity function wrt F"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "psi = strainDensityFunctionNeoHookean(F=F, Ey=1000, nu=0.4)\n",
+ "\n",
+ "integral_E = psi * fe.dx\n",
+ "\n",
+ "Pi = integral_E - integral_S\n",
+ "\n",
+ "# Compute 1st variation of Pi (directional derivative about u in dir. of v)\n",
+ "Fpi = fe.derivative(Pi, u, v)\n",
+ "\n",
+ "# Compute Jacobian of F\n",
+ "J = fe.derivative(Fpi, u, du)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(4, True)"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "############################## SOLVER PARAMS ##################################\n",
+ "# Define the solver params\n",
+ "problem = fe.NonlinearVariationalProblem(Fpi, u, bcs, J)\n",
+ "solver = fe.NonlinearVariationalSolver(problem)\n",
+ "prm = solver.parameters\n",
+ "prm['nonlinear_solver'] = 'newton'\n",
+ "prm['newton_solver']['linear_solver'] = 'petsc'\n",
+ "\n",
+ "prm['newton_solver']['error_on_nonconvergence'] = True\n",
+ "prm['newton_solver']['absolute_tolerance'] = 1E-9\n",
+ "prm['newton_solver']['relative_tolerance'] = 1E-8\n",
+ "prm['newton_solver']['maximum_iterations'] = 25\n",
+ "prm['newton_solver']['relaxation_parameter'] = 1.0\n",
+ "\n",
+ "prm['newton_solver']['lu_solver']['report'] = True\n",
+ "#prm['newton_solver']['lu_solver']['reuse_factorization'] = False\n",
+ "#prm['newton_solver']['lu_solver']['same_nonzero_pattern'] = False\n",
+ "prm['newton_solver']['lu_solver']['symmetric'] = False\n",
+ "\n",
+ "prm['newton_solver']['krylov_solver']['error_on_nonconvergence'] = True\n",
+ "prm['newton_solver']['krylov_solver']['absolute_tolerance'] = 1E-7\n",
+ "prm['newton_solver']['krylov_solver']['relative_tolerance'] = 1E-5\n",
+ "prm['newton_solver']['krylov_solver']['maximum_iterations'] = 1000\n",
+ "prm['newton_solver']['krylov_solver']['nonzero_initial_guess'] = True\n",
+ "if prm['newton_solver']['linear_solver'] == 'gmres':\n",
+ " prm['newton_solver']['preconditioner'] = 'ilu'\n",
+ "\n",
+ "solver.solve()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "boundary_file = fe.File(\"result_growth.pvd\")\n",
+ "boundary_file << u"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Saving displacement solution to file...\n"
+ ]
+ }
+ ],
+ "source": [
+ "\n",
+ "\n",
+ "############################# POST-PROCESSING ################################\n",
+ "# Save solution to file in VTK format\n",
+ "print(\"Saving displacement solution to file...\")\n",
+ "uViewer = fe.File('paraview/displacement.pvd')\n",
+ "uViewer << u\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Min/Max displacement: 5.255531433076261e-06 1.0846528324258178\n"
+ ]
+ }
+ ],
+ "source": [
+ "W = fe.FunctionSpace(mesh, 'P', 1)\n",
+ "\n",
+ "# Maximum and minimum displacement\n",
+ "u_magnitude = fe.sqrt(fe.dot(u, u))\n",
+ "u_magnitude = fe.project(u_magnitude, W)\n",
+ "print('Min/Max displacement:',\n",
+ " u_magnitude.vector().get_local().min(),\n",
+ " u_magnitude.vector().get_local().max())\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Computing the deformation tensor and saving to file...\n"
+ ]
+ }
+ ],
+ "source": [
+ "Z = fe.TensorFunctionSpace(mesh, 'P', 1)\n",
+ "\n",
+ "# Computation of the large deformation strains\n",
+ "print(\"Computing the deformation tensor and saving to file...\")\n",
+ "epsilon_u = largeKinematics(u)\n",
+ "epsilon_u_project = fe.project(epsilon_u, Z)\n",
+ "epsilonViewer = fe.File('paraview/strain.pvd')\n",
+ "epsilonViewer << epsilon_u_project\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Stress derivation and saving to file...\n"
+ ]
+ }
+ ],
+ "source": [
+ "\n",
+ "# Computation of the stresses\n",
+ "print(\"Stress derivation and saving to file...\")\n",
+ "S = fe.diff(psi, F)\n",
+ "S_project = fe.project(S, Z)\n",
+ "sigmaViewer = fe.File('paraview/stress.pvd')\n",
+ "sigmaViewer << S_project"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Maximum equivalent stress: 118.5307998241279\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Computation of an equivalent stress\n",
+ "s = S - (1./3)*fe.tr(S)*fe.Identity(u.geometric_dimension())\n",
+ "von_Mises = fe.sqrt(3./2*fe.inner(s, s))\n",
+ "von_Mises_project = fe.project(von_Mises, W)\n",
+ "misesViewer = fe.File('paraview/mises.pvd')\n",
+ "misesViewer << von_Mises_project\n",
+ "print(\"Maximum equivalent stress:\", von_Mises_project.vector().get_local().max())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/FEMII/project_growth_model_fenics.ipynb b/FEMII/project_growth_model_fenics.ipynb
index c83852c..4eb7ae9 100644
--- a/FEMII/project_growth_model_fenics.ipynb
+++ b/FEMII/project_growth_model_fenics.ipynb
@@ -1,94 +1,72 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
- "# Finite Growth \n",
- "using multiplicative Decomposition of Deformation Gradient\n",
+ "# Finite Growth (FEniCS)\n",
"\n",
"*Rafael Gfeller*\n",
"\n",
- "# Preparation:\n",
- "- theory overview\n",
- "- simple example formulation for F_g\n",
+ "## Tasks\n",
"\n",
+ "The goal of this project is to implement and evaluate a finite growth mechanism based on the multiplicative decomposition of the \n",
+ "deformation gradient.\n",
+ "The model will be implemented in FEniCS and the implementation will be verified by comparison to literature results.\n",
+ "We then compare results of this (non-linear) growth model to a linearized version, both in test geometries and in\n",
+ "a realistic model of human brain anatomy.\n",
"\n",
+ "- familiarize with FEniCS: [Notebook & Links](https://github.com/danielabler/femII_fenics_2019/blob/master/00_intro_to_fenics.ipynb) \n",
+ "- familiarize with hyperelasticity and implementation in FEniCS \n",
+ "- familiarize with finite growth theory\n",
+ "- implement finite growth model in [Feng (2018)](https://doi.org/10.1007/s00466-018-1589-2), Appendix A, and linearized growth model, equations (1)-(7)\n",
+ "- compare both growth models, [Feng (2018)](https://doi.org/10.1007/s00466-018-1589-2), Fig. 13, 14\n",
+ "- apply model to brain anatomy \n",
"\n",
+ "Additional Links:\n",
+ "- [Feng (2018)](https://doi.org/10.1007/s00466-018-1589-2)\n",
"\n",
- "s you have seen in class yesterday, the goal is to develop, implement and \n",
- "test a finite growth mechanism based on the multiplicative decomposition of the \n",
- "deformation gradient. The idea behind this approach is to decompose \n",
- "deformations taking place due to a growth process into a growth-related and an \n",
- "elastic component. The elastic part is determined by the material model, the \n",
- "growth related part is driven by some internal (and problem specific) process. \n",
- "I am interested in the application of this model for tumor growth simulation, \n",
- "so the growth component could for example be described as a function of cell \n",
- "density.\n",
- "\n",
- "You would first be familiarizing with FENICS (https://fenicsproject.org/) , for \n",
- "example by implementing a simple hyperelastic material model such as this \n",
- "tutorial example: \n",
- "https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/\n",
- "hyperelasticity/python/documentation.html .\n",
- "\n",
- "You would then extend this example to include a simple growth mechanism. To \n",
- "implement this you will need to learn and understand the theory behind. \n",
- "\n",
- "Once this works, we can decide whether to investigate different growth \n",
- "mechanisms (for tumor growth, possibly for cornea with Miguel), or to combine \n",
- "this mechanism with different hyperelastic material models.\n",
- "\n",
- "You will initially be working on very simple test geometries. If progress \n",
- "allows, you can apply the model in the end to more realistic geometries, such \n",
- "as a model of the human brain.\n",
- "\n",
- "FENICS is a finite element modeling framework with a python interface. A nice \n",
- "feature (which I find very instructive) of this framework is that the \n",
- "implementation is very close to the mathematics of the underlying equations. \n",
- "There is no GUI, so yes, some programming will be involved for implementation, \n",
- "as well as for post processing of simulation results. Everything can be done \n",
- "in python.\n",
- "\n",
- "Let me know if you have more questions or would like any other information.\n",
- "\n",
- "best,\n",
- "Daniel"
+ "## You will learn\n",
+ "- FEniCS & Python\n",
+ "- hyperelasticity theory\n",
+ "- finite growth theory\n",
+ "- plane strain / plane stress\n",
+ "- application to tumor growth modeling\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"pycharm": {
"stem_cell": {
"cell_type": "raw",
"source": [],
"metadata": {
"collapsed": false
}
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
\ No newline at end of file
diff --git a/FEMII/project_hyperelastic_model_abaqus.ipynb b/FEMII/project_hyperelastic_model_abaqus.ipynb
index 082fb94..997a653 100644
--- a/FEMII/project_hyperelastic_model_abaqus.ipynb
+++ b/FEMII/project_hyperelastic_model_abaqus.ipynb
@@ -1,60 +1,70 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"\n",
- "# Hyperelastic Material Model for Brain Tissue\n",
+ "# Hyperelastic Material Model for Brain Tissue (Abaqus)\n",
"\n",
"*Marc Sascha Ilic*\n",
"\n",
"\n",
"## Tasks:\n",
"- implement Abaqus Ogden model \n",
"- reproduce stress-shear/stretchrelations from Fig 24 in in simple geometry, using parameter values from Table 1 in [1]\n",
"- optimize using FEM model -- better fit possible?\n",
"- linearize and compare linearized and non linear models\n",
"- calibrate linear elastic model to best agree with experimental stress/strain relation in compression\n",
"- enable thermal expansion in linear elastic and hyperelastic material; compare temporal stress evolution in test geometry\n",
"- apply both models to few seeding scenarios in brain anatomy and compare resulting tumor shapes and stress distribution.\n",
"\n",
- "[1] https://link.springer.com/content/pdf/10.1007%2Fs11831-019-09352-w.pdf \n"
+ "\n",
+ "Additional Links:\n",
+ "- [Ogden (1971)](https://www.jstor.org/stable/77930)\n",
+ "- [Budday (2017)](https://www.sciencedirect.com/science/article/pii/S1742706116305633)\n",
+ "- [Budday (2019)](https://link.springer.com/article/10.1007%2Fs11831-019-09352-w)\n",
+ "\n",
+ "## You will learn\n",
+ "- Abaqus\n",
+ "- hyperelasticity theory\n",
+ "- characterization of material models\n",
+ "- material model parameter optimization"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"pycharm": {
"stem_cell": {
"cell_type": "raw",
"source": [],
"metadata": {
"collapsed": false
}
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
\ No newline at end of file
diff --git a/FEMII/project_hyperelastic_model_fenics.ipynb b/FEMII/project_hyperelastic_model_fenics.ipynb
index d45fc5e..ae1f0fb 100644
--- a/FEMII/project_hyperelastic_model_fenics.ipynb
+++ b/FEMII/project_hyperelastic_model_fenics.ipynb
@@ -1,64 +1,72 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
- "# Hyperelastic Material Model for Brain Tissue\n",
+ "# Hyperelastic Material Model for Brain Tissue (FEniCS)\n",
"\n",
"*Mathieu Simon*\n",
"\n",
"## Tasks:\n",
"\n",
"The goal of this project is to implement and evaluate a hyperelastic (Ogden) material model for brain tumor growth simulations.\n",
"The model will be implemented in FEniCS and the implementation will be verified by comparison to results of conventional FEM software, such as Abaqus,\n",
"and validated by comparison to published experimental results.\n",
"We will then use this model to evaluate tissue deformations and stresses caused by [hydrocephalus](https://en.wikipedia.org/wiki/Hydrocephalus) \n",
"in a realistic model of human brain anatomy, and compare to results using a linear elastic material model.\n",
"\n",
- "- familiarize with FEniCS\n",
- "- implement Ogden model, eq. 8 in [1] \n",
- "- reproduce stress-shear/stretchrelations from Fig 24 in simple geometry, using parameter values from Table 1 in [1]\n",
+ "- familiarize with FEniCS: [Notebook & Links](https://github.com/danielabler/femII_fenics_2019/blob/master/00_intro_to_fenics.ipynb) \n",
+ "- familiarize with hyperelasticity and implementation in FEniCS \n",
+ "- implement Ogden model, eq. 8 in [Budday (2019)](https://link.springer.com/article/10.1007%2Fs11831-019-09352-w)\n",
+ "- using parameter values from Table 1 in [Budday (2019)](https://link.springer.com/article/10.1007%2Fs11831-019-09352-w), reproduce stress-shear/stretchrelations from Fig 24\n",
"- apply model to brain anatomy \n",
"\n",
- "- Ogden, 1971: https://www.jstor.org/stable/77930\n",
- "- Budday 2017: https://www.sciencedirect.com/science/article/pii/S1742706116305633\n",
- "- [1] Budday 2019: https://link.springer.com/article/10.1007%2Fs11831-019-09352-w"
+ "Additional Links:\n",
+ "- [Ogden (1971)](https://www.jstor.org/stable/77930)\n",
+ "- [Budday (2017)](https://www.sciencedirect.com/science/article/pii/S1742706116305633)\n",
+ "- [Budday (2019)](https://link.springer.com/article/10.1007%2Fs11831-019-09352-w)\n",
+ "\n",
+ "## You will learn\n",
+ "- FEniCS & Python\n",
+ "- hyperelasticity theory\n",
+ "- characterization of material models\n",
+ "- brain tissue modeling "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"pycharm": {
"stem_cell": {
"cell_type": "raw",
"source": [],
"metadata": {
"collapsed": false
}
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
\ No newline at end of file