nonlinear_heat_transfer.ipynb 36.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stationnary non-linear heat transfer\n",
    "\n",
    "## Description of the non-linear constitutive heat transfer law\n",
    "\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
    "The thermal material is described by the following non linear Fourier\n",
    "Law:\n",
    "\n",
    "$$\n",
    "\\mathbf{j}=-k\\left(T\\right)\\,\\mathbf{\\nabla} T\n",
    "$$\n",
    "\n",
    "where $\\mathbf{j}$ is the heat flux and $\\mathbf{\\nabla} T$ is the\n",
    "temperature gradient.\n",
    "\n",
    "Expression of the thermal conductivity\n",
    "--------------------------------------\n",
    "\n",
    "The thermal conductivity is assumed to be given by:\n",
    "\n",
    "$$\n",
    "k\\left(T\\right)={\\displaystyle \\frac{\\displaystyle 1}{\\displaystyle A+B\\,T}}\n",
    "$$\n",
    "\n",
    "This expression accounts for the phononic contribution to the thermal\n",
    "conductivity.\n",
    "\n",
    "Derivatives\n",
    "-----------\n",
    "\n",
    "As discussed below, the consistent linearisation of the heat transfer\n",
    "equilibrium requires to compute:\n",
    "\n",
    "-   the derivative\n",
    "    ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}$\n",
    "    of the heat flux with respect to the temperature gradient.\n",
    "    ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}$\n",
    "    is given by: $$\n",
    "      {\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}=-k\\left(T\\right)\\,\\matrix{I}\n",
    "    $$\n",
    "-   the derivative\n",
    "    ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}$\n",
    "    of the heat flux with respect to the temperature.\n",
    "    ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}$\n",
    "    is given by: $$\n",
    "      {\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}=-{\\displaystyle \\frac{\\displaystyle \\partial k\\left(T\\right)}{\\displaystyle \\partial T}}\\,\\mathbf{\\nabla} T=B\\,k^{2}\\,\\mathbf{\\nabla} T\n",
    "    $$\n",
    "\n",
    "`MFront`’ implementation\n",
    "========================\n",
    "\n",
    "Choice of the the domain specific language\n",
    "------------------------------------------\n",
    "\n",
    "Every `MFront` file is handled by a domain specific language (DSL), which\n",
    "aims at providing the most suitable abstraction for a particular choice\n",
    "of behaviour and integration algorithm. See `mfront mfront --list-dsl`\n",
    "for a list of the available DSLs.\n",
    "\n",
    "The name of DSL’s handling generic behaviours ends with\n",
    "`GenericBehaviour`. The first part of a DSL’s name is related to the\n",
    "integration algorithm used.\n",
    "\n",
    "In the case of this non linear transfer behaviour, the heat flux is\n",
    "explicitly computed from the temperature and the temperature gradient.\n",
    "The `DefaultGenericBehaviour` is the most suitable choice:\n",
    "\n",
    "``` cxx\n",
    "@DSL DefaultGenericBehaviour;\n",
    "```\n",
    "\n",
    "Some metadata\n",
    "-------------\n",
    "\n",
    "The following lines define the name of the behaviour, the name of the\n",
    "author and the date of its writing:\n",
    "\n",
    "``` cxx\n",
    "@Behaviour StationaryHeatTransfer;\n",
    "@Author Thomas Helfer;\n",
    "@Date 15/02/2019;\n",
    "```\n",
    "\n",
    "Gradients and fluxes\n",
    "--------------------\n",
    "\n",
    "Generic behaviours relate pairs of gradients and fluxes. Gradients and\n",
    "fluxes are declared independently but the first declared gradient is\n",
    "assumed to be conjugated with the first declared fluxes and so on…\n",
    "\n",
    "The temperature gradient is declared as follows (note that Unicode characters are supported):\n",
    "\n",
    "``` cxx\n",
    "@Gradient TemperatureGradient ∇T;\n",
    "∇T.setGlossaryName(\"TemperatureGradient\");\n",
    "```\n",
    "\n",
    "Note that we associated to `∇T` the glossary name `TemperatureGradient`.\n",
    "This is helpful for the calling code.\n",
    "\n",
    "After this declaration, the following variables will be defined:\n",
    "\n",
    "-   The temperature gradient `∇T` at the beginning of the time step.\n",
    "-   The increment of the temperature gradient `Δ∇T` over the time step.\n",
    "\n",
    "The heat flux is then declared as follows:\n",
    "\n",
    "``` cxx\n",
    "@Flux HeatFlux j;\n",
    "j.setGlossaryName(\"HeatFlux\");\n",
    "```\n",
    "\n",
    "In the following code blocks, `j` will be the heat flux at the end of\n",
    "the time step.\n",
    "\n",
    "Tangent operator blocks\n",
    "-----------------------\n",
    "\n",
    "By default, the derivatives of the gradients with respect to the fluxes\n",
    "are declared. Thus the variable `∂j∕∂Δ∇T` is automatically declared.\n",
    "\n",
    "However, as discussed in the next section, the consistent linearisation\n",
    "of the thermal equilibrium requires to return the derivate of the heat\n",
    "flux with respect to the increment of the temperature (or equivalently\n",
    "with respect to the temperature at the end of the time step).\n",
    "\n",
    "``` cxx\n",
    "@AdditionalTangentOperatorBlock ∂j∕∂ΔT;\n",
    "```\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "\n",
    "The `A` and `B` coefficients that appears in the definition of the\n",
    "thermal conductivity are declared as parameters:\n",
    "\n",
    "``` cxx\n",
    "@Parameter real A = 0.0375;\n",
    "@Parameter real B = 2.165e-4;\n",
    "```\n",
    "\n",
    "Parameters are stored globally and can be modified from the calling\n",
    "solver or from `python` in the case of the coupling with `FEniCS`\n",
    "discussed below.\n",
    "\n",
    "Local variable\n",
    "--------------\n",
    "\n",
    "A local variable is accessible in each code blocks.\n",
    "\n",
    "Here, we declare the thermal conductivity `k` as a local variable in\n",
    "order to be able to compute its value during the behaviour integration\n",
    "and to reuse this value when computing the tangent operator.\n",
    "\n",
    "``` cxx\n",
    "@LocalVariable thermalconductivity k;\n",
    "```\n",
    "\n",
    "Integration of the behaviour\n",
    "----------------------------\n",
    "\n",
    "The behaviour integration is straightforward: one starts to compute the\n",
    "temperature at the end of the time step, then we compute the thermal\n",
    "conductivity (at the end of the time step) and the heat flux using the\n",
    "temperature gradient (at the end of the time step).\n",
    "\n",
    "``` cxx\n",
    "@Integrator{\n",
    "  // temperature at the end of the time step\n",
    "  const auto T_ = T + ΔT;\n",
    "  // thermal conductivity\n",
    "  k = 1 / (A + B ⋅ T_);\n",
    "  // heat flux\n",
    "  j = -k ⋅ (∇T + Δ∇T);\n",
    "} // end of @Integrator\n",
    "```\n",
    "\n",
    "Tangent operator\n",
    "----------------\n",
    "\n",
    "The computation of the tangent operator blocks is equally simple:\n",
    "\n",
    "``` cxx\n",
    "@TangentOperator {\n",
    "  ∂j∕∂Δ∇T = -k ⋅ tmatrix<N, N, real>::Id();\n",
    "  ∂j∕∂ΔT  =  B ⋅ k ⋅ k ⋅ (∇T + Δ∇T);\n",
    "} // end of @TangentOperator \n",
    "```\n",
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396
    "## FEniCS implementation\n",
    "\n",
    "We consider a rectanglar domain with imposed temperatures `Tl` (resp. `Tr`) on the left (resp. right boundaries). We want to solve for the temperature field `T` inside the domain using a $P^1$-interpolation. We initialize the temperature at value `Tl` throughout the domain. We finally load the material library with a `plane_strain` hypothesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from dolfin import *\n",
    "import mfront_wrapper as mf\n",
    "\n",
    "length = 30e-3\n",
    "width = 5.4e-3\n",
    "mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)\n",
    "\n",
    "V = FunctionSpace(mesh, \"CG\", 1)\n",
    "T = Function(V, name=\"Temperature\")\n",
    "\n",
    "def left(x, on_boundary):\n",
    "    return near(x[0], 0) and on_boundary\n",
    "def right(x, on_boundary):\n",
    "    return near(x[0], length) and on_boundary\n",
    "\n",
    "Tl = 300\n",
    "Tr = 800\n",
    "T.interpolate(Constant(Tl))\n",
    "\n",
    "bc = [DirichletBC(V, Constant(Tl), left),\n",
    "      DirichletBC(V, Constant(Tr), right)]\n",
    "\n",
    "material = mf.MFrontNonlinearMaterial(\"../materials/src/libBehaviour.so\",\n",
    "                                      \"StationaryHeatTransfer\",\n",
    "                                      hypothesis=\"plane_strain\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MFront behaviour implicitly declares the temperature as an external state variable called `\"Temperature\"`. We must therefore associate this external state variable to a known mechanical field. This can be achieved explicitly using the `register_external_state_variable` method. In the present case, this can be donc automatically since the name of the unknown temperature field matches the [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) name `\"Temperature\"`. In this case, the following message is printed:\n",
    "```\n",
    "Automatic registration of 'Temperature' as an external state variable.\n",
    "```\n",
    "For problems in which the temperature only acts as a parameter (no jacobian blocks with respect to the temperature), the temperature can be automatically registered as a constant value ($293.15 \\text{ K}$ by default) or to any other (`dolfin.Constant` or `float`) value using the `register_external_state_variable` method.\n",
    "\n",
    "In the FEniCS interface, we instantiate the main mechanical unknown, here the temperature field `T` which has to be named `\"Temperature\"` in order to match MFront's predefined name. Using another name than this will later result in an error saying:\n",
    "```\n",
    "ValueError: 'Temperature' could not be associated with a registered gradient or a known state variable.\n",
    "```\n",
    "\n",
    "The MFront behaviour declares the field `\"TemperatureGradient\"` as a Gradient variable, with its associated Flux called `\"HeatFlux\"`. We can check that the `material` object retrieves MFront's gradient and flux names, as well as the different tangent operator blocks which have been defined, namely `dj_ddgT` and `dj_ddT` in the present case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['TemperatureGradient']\n",
      "['HeatFlux']\n",
      "['dHeatFlux_dTemperatureGradient', 'dHeatFlux_dTemperature']\n"
     ]
    }
   ],
   "source": [
    "print(material.get_gradient_names())\n",
    "print(material.get_flux_names())\n",
    "print([\"d{}_d{}\".format(*t) for t in material.get_tangent_block_names()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When defining the non-linear problem, we will specify the boundary conditions and the requested quadrature degree which will control the number of quadrature points used in each cell to compute the non-linear constitutive law. Here, we specify a quadrature of degree 2 (i.e. 3 Gauss points for a triangular element). Finally, we need to associate to MFront gradient object the corresponding UFL expression as a function of the unknown field `T`. To do so, we use the `register_gradient` method linking MFront `\"TemperatureGradient\"` object to the UFL expression `grad(T)`. Doing so, the corresponding non-linear variational problem will be automatically be built:\n",
    "\n",
    "\\begin{equation}\n",
    "F(\\widehat{T}) = \\int_\\Omega \\boldsymbol{j}\\cdot \\nabla \\widehat{T} \\text{dx} = 0 \\quad \\forall \\widehat{T}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=2, bcs=bc)\n",
    "problem.register_gradient(\"TemperatureGradient\", grad(T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the two tangent operator blocks `dj_ddgT` and `dj_ddT`, it will automatically be deduced that the heat flux $\\boldsymbol{j}$ is a function of both the temperature gradient $\\boldsymbol{g}=\\nabla T$ and the temperature itself i.e. $\\boldsymbol{j}=\\boldsymbol{j}(\\boldsymbol{g}, T)$. The following tangent bilinear form will therefore be used when solving the above non-linear problem:\n",
    "\n",
    "\\begin{equation}\n",
    "J(\\widehat{T},T^*) = \\int_{\\Omega} \\nabla \\widehat{T}\\cdot\\left(\\dfrac{\\partial \\boldsymbol{j}}{\\partial \\boldsymbol{g}}\\cdot \\nabla T^*+\\dfrac{\\partial \\boldsymbol{j}}{\\partial T}\\cdot T^*\\right) \\text{dx}\n",
    "\\end{equation}\n",
    "\n",
    "Similarly to the case of external state variables, common gradient expressions for some [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) names have been already predefined which avoid calling explicitly the `register_gradient` method. Predefined expressions can be obtained from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'TFEL gradient name'   (Available hypotheses)\n",
      "---------------------------------------------\n",
      "'Strain'               ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')\n",
      "'TemperatureGradient'  ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')\n",
      "'DeformationGradient'  ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mf.list_predefined_gradients()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the name `\"Temperature Gradient\"` is in fact a predefined gradient. Omitting calling the `register_gradient` method will in this case print the following message upon calling `solve`:\n",
    "```\n",
    "Automatic registration of 'TemperatureGradient' as grad(Temperature).\n",
    "```\n",
    "meaning that a predefined gradient name has been found and registered as the UFL expression $\\nabla T$.\n",
    "\n",
    "We finally solve the non-linear problem using a default Newton non-linear solver. The `solve` method returns the number of Newton iterations (4 in the present case) and converged status ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatic registration of 'Temperature' as an external state variable.\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(4, True)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "problem.solve(T.vector())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finally check that the thermal conductivity coefficient $k$, computed from the ratio between the horizontal heat flux and temperature gradient matches the temperature-dependent expressions implemented in the MFront behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEUCAYAAAA8+dFZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvW9wVGeWp/kcSSALBJZEDdAFRAthTBdQIVuSW0sMZfwHV82GIzZ63eCadWzEdMwa8HyeKoN7NmK/7CyWqz5vD1CxMfOlOmzjmt6NdURXG7vBzQRDW5LNGKgGjMQGZhqYQsIgLGQknf1w7833KpX/82bem5nnichQ5s3Me49+efOefN/znnNEVTEMwzCMQmmK2wDDMAyjtjDHYRiGYRSFOQ7DMAyjKMxxGIZhGEVhjsMwDMMoCnMchmEYRlGY4zAMwzCKouqOQ0Q+yrCtR0TeFJHd/t+OattlGIZhFIZUKwFQRHYDPcARVZW050ZUtd+/3wEcU9W9VTHMMAzDKIqqOY7UAUUmVbUz9LgPGFLVl7K9xjAMw0gOSYhx9AB307ZNiEhPHMYYhmEYuWmJ2wCgK8v2jHEOEdkP7Adoa2vr37BhAwCtra00Nzfz7bffAtDS0kJbWxv3798P3kd7ezvffvstc3NzACxfvpxHjx7x3XffAfDYY48hIkxPTwOwZMkSWltbmZqaAqCpqYnly5dHso8HDx4wPz8PQHt7OzMzMzx69Aj//0JVefjwIQBLly5lyZIlPHjwAIDm5maWLVu2YB9Lly5FVYvex9TUFMGoc8WKFUxPTzM7OwvAsmXLmJubY2ZmpiCNo9hHFJ9TYEcSP6dS91Hq56SqLFmyJJGfU7W/T9PT04hIIj+nTBpX6nO6fPny71X1n1AGSXAcEyx2EtmcCap6FDgKsGXLFr106VIFTasdTp48yXPPPRe3GYnAtHCYFg7TwkNE/r9y91GU4xCRbmAv8BIQDo7cBd5V1d+UYMNYpo2qOprvjcuWLSvhcPVJb29v3CYkBtPCYVo4TIvoKCjGISKPi8jbwG7guKr+WFV/Err9FPhGRP6diLxSjAG+g0iNMPzYxolC3hsMvQwWDG8bHdPCYVo4TIvoyOs4RORxYK+qHlLVX6nqeKbXqerHqvoG8Hkm5yEifSLyJtAhIkP+8tyAfUEeB7AH2FeI8cE8nwFjYxkHbg2JaeEwLRymRXTknapS1W+AXxW6Q9+xLHIu/shiFHgnx3NQ4GjDMAzDiIckLMctmWv35uk+9CHdhz6M25TY6e7ujtuExGBaOEwLh2kRHSUlAIrIC6r6SQXsKYqB7zfrf/6X7Twx+2sArr39cswWxce9e/dYuXJl3GYkAtPC0Yha3Lt3j9u3b6eW0gbMz8/T1FTTv5ULYvny5axfvz7r/+pX6hgo5xg5p6r8+MYIcDW8GdgIbC7nwFHR3Axf8RpPzP6a7kMfNqzzGB0dtaWGPqaFo9G0uHfvHrdu3WLdunW0tbUh4qob3b9/nxUrVsRoXeWZn5/nxo0b/P73v2f16tUVO05O9+vHN/amraD6MfBGxSwqgaYm+KrlNQCbtjKMBub27dusW7eOZcuWLXAajUJTUxNr1qzhm2++qexx8r1AVT/PsO3jyphTGiL+yKOBnUdnp5X2CjAtHI2mxaNHj2hra8v4XHNzc5WtiYclS5akMs0rRc1P+Il4N2hs52HJTQ7TwtGIWmQbaTRKwnA1RlpFOY5ik/uqQRDbDzuP4ZaC0kDqilOnTsVtQmIwLRymhcMSAKOj2BFH1hpScRA4DVXvFjiPriav8FgjjTqqXR4/yZgWDtMi+bz00kuICJ2dnQtuQcJituePHz8em801PVV1nh7m5pwDCRBpvCmrRgwEZsO0cJgWtcHQ0BCTk5MLbj09PTmf37NnT2z21rTjAHhi9tf4lZAXjDoaLd6xa9euuE1IDKaFw7Rw1PtS3GpS047jia6l3t/ZXxPUO2xU53Hu3Lm4TUgMpoXDtHAEPSsK4a8+v8E/ffsTNh76kH/69if81ec3KmhZ7VGs40jUuHd2djaV8NfozmNycjJuExKDaeEwLRyFVtP+q89v8NZvvuTG3WkUuHF3mrd+86U5jxBFOQ5VPVYpQ8rBnIdhGFHxi99eYvrRQicz/WiOX/y2ck3jDh48uCDwffDgwZzPx52fU7DjEJHDlTSkFMLrshvdefT19cVtQmIwLRymhaPQPI7/ene6qO1RkB78Hhoayvl83CPJYkYcBypmRYmkDz0b2XlMTEzEbUJiMC0cpoWj0Gzq73dkzjzPtr0RKcZx7M2VACgiP4vAnqLI1MipUZ3HtWvX4jYhMZgWDtPC8d133xX0up//ZAttSxaWJ2lb0szPf7KlEmbVJAU7Dr8+1cci8lR4u4g8JSLvAUOZ31l9GtV5GIZRPn/y9DoOv/JD1nW0IcC6jjYOv/JD/uTpdXGblhiKDY5/A4iIdIvIKyIyDHwMfAZUPXDe2tqa9blGcx7hZKFGx7RwmBaOpUuXFvzaP3l6Hf/p0AuMv/0y/+nQC7E7jYMHDyIiC27pAfRqUnAjJxF5SlW/8O+/Dfwp8E6cK62efvpp/fzzRcV7FxA4hbGlry0oiBj823Nz1EUjqMnJydhXWiQF08LRaFr87ne/4wc/+EHG52ZnZ2lpydstuy7IpUMUjZxKCo6r6iHgjbiX5xaS0BM4g3B2eT2OPCzRy2FaOEwLx/R05VZFNRpFOQ4RmROR3/qB8MfD8Q4RSXR/ymDKKjzAqjfnYRiGUQ2KHXF0AUeBJ4BfASMi8q6IvE4MwfFCh53heMfE/HKg/kYeq1atituExGBaOEwLR6M0cqoGxayqOqaq36jqB6r6hqp24TmQE8CrwP5KGZmNbJ2+MhE4j4HZY9yZW76gJHs9OI9t27bFbUJiMC0cpoWjmOuFkZuyihyq6rjvUH4MHIrIpoIptjFL2Hn0fFdfq60+/fTTuE1IDKaFw7RwTE1NxW1C3RBlddyjEe6rYlx7++UFU1cW8zAMwyiOyByHn+NRVcppUhM4j3qJeTTKMsNCMC0cpoVRCSJxHCLyuB8gryrt7e1l7yM95gG16Tx27twZtwmJwbRwmBYOa+QUHVGNOLqAH5ezAxHpE5E3RWSP/7cj33uKacySiXDMo9bzPEZHR+M2ITGYFg7TwvHgwYO4TagbInEcfpD81VLf7zuJ91X1HVU9DhyngOW9hTZmyUV6aZJaXW117969uE1IDKaFw7RwzAe/DmucAwcO0NnZyd69e2OzISmtY3cDY8EDVR2jist7w86j3lZbGYZRPxw96q1Bmpyc5Nix+Ap3lOQ4ROSFiO24izfdlX6cnBXali9fHpkBtb7aamCgrNIzdYVp4TAtHIU2ckoyIyMj9Pf3A9DRkXc2v2LkXHIhIo8DI8DV8GZgI7A5KiNU9YRf8bFDVe+KSNC2bJEyIrIffzSydu1aTp48CXhVQFesWJGqzbNq1Sq2bduWWsfe0tLCzp07GR0dTQ3fBwYGuHXrFtevXwdg8+bNDP9skIFfnmV+3nMUgQMR8e43N8NXvMYTs7+m+9CH/O5/e4Gvv/6aGze8fsRbtmyhubmZixcvEti4ceNGzpw5A3hJSIODg5w9ezZVO2fHjh2Mj49z8+ZNALZu3crc3ByXLnmtKtetW8f69es5e/Ys4C0KGBgY4MyZM6meJN///veZnZ3l9u3bAGzfvp2ZmRmuXLkCwIYNG1izZg3Dw8MArFy5kr6+Pk6fPp1qcPPss89y4cIF7ty5A0Bvby/3799nbMwbDHZ3d9PV1ZWaN+/s7KS3t5dTp06hqogIu3bt4ty5c6kOZX19fUxMTKT6QkT1ObW2tnL+/HkAVq9ezZNPPsnp06cBbwrzxRdfZHh4OLV2f3BwMBGf086dO7l8+XLVPqeZmRnWrl2byM+ptbWVHTt2RPo5Pf7448zPzzMzM5PS67HHHkNVmZ6epqmpiSVLlrB06dJUzKOpqYnly5czNTVFUPS1vb2dhw8fLtpH8DkuWbKEJUuWpOKswT7CuWXt7e1MT0+nptTb2tqYm5tL9QVZunQpLS0tqX00NzezbNmyBftYsWIF3377bWofqsqjR49Sr1m6dCnNzc2pc7S5uZm2tjYePnzIyZMnM35OUZC3Oq6IPK2qn6dte9HvzxEpIrIHb8rqLp6z6lTVu9lev2XLFg2+tFESjCa+anmNpiY34gicByysqgvxV9Y9efIkzz33XKw2JAXTwtFoWuSqCnv//v3iVlZd/3u49nfQ/SPY8McRWbiYzs5OPv74Y/bt28eBAwfYv38/o6Oj7Nu3j7GxMfbv38/Q0BAHDx7knXfeoaOjg56eHkZGRrLuM/bquOlOw98WudPw93tcVYNlIGO5nEYlKTTmMbb0NYZb9gHJm7oyDKMMrv89/If/AT75t97f639f0cPt3buXoaEh9u/3Qrv79u3j/fffT40Ijx49mnp+aGgop9OoBkkJjiMi4e7rB4C8XUoee+yxitmTHvPI5DxEYFXzg0TEPTZvjmzmsOYxLRymhSNX47dFXPs7mPsOdM77e+3vKmcYXqOm3bt3A56TGBsbY+/evfT393P8+PHYHUU6RTmOXD3HI+Cgn8OxH/jMX5abz54KmuORy3kkKWhe1JeizjEtHKaFo6jrRfePoHkpSLP3t/tHlTMMUk4D4OrVq7z11luMjIwwMjLC1atXOXLkSEWPXyzFjjgWrXyKClU96k9VHS3EaUD1GrNkcx7ZluvGQRCANEyLMKaF4+HDh4W/eMMfw7/4f+CFf+P9rWCMA6Cry11an3nmGT766KOKHq9cEjNVlXSyJQqCcx5Nvprdhz60mIdh1Dob/hh+9K8r7jTS2bNnD8PDwxw/7v1+Hh0dTd1PCjXtOJYsWVLV44WdR7hESYBIfAHz1atXV+1YSce0cJgWjloq+DgyMsLhw4fp7Oxk37599PTkTGmrOnmX4y54scjrqvqrCtpTFP39/RpH0Ci8XDdoKhZeqgsLl+tWY6nu7OxsTX0xKolp4Wg0LXItQw1yWBqB2Jfjph+znINFTVyNWUoJmFd69BEkVRmmRRjTwmGNnKKjKMehqvEVR0kYhQbMx5a+FvuKK8MwjCiJqh/Hi1Hsp1iamuIN0RRSWVekOst1bdmlw7RwmBaORpmmqgalFjncKCI/CxU7HK5wjkdGoixyWCr5sszDK64qGTTfsWNH5PusVUwLh2nhiKLxm+FR6k/2PcA48KqIXMHrN/5MZFYVSLmNnKIiV5Z5QKWzzIOieIZpEca0cFgjp+go1XGMquoHqvqGqm4G3sZzHlUlikZOURJnlrkF/hymhaMRtci2UrReGjnlo5iVsqVScpBARFYG91X1c1Udj8ak2ibXyCM8bWUBc8OIniVLllStokRSefToUcWXYJcTXR4VkcMVaOpUMEmIcWSikGmrqEceg4ODZe+jXjAtHI2mxerVq7lx4wbffvvtol/eSb1eRMn8/Dy3bt3i8ccfr+hxikoATL1JZB9wAugDXsJr/Tqiqj+N1rzc/PCHP9Qvv/yymocsimyJghB9X48rV65YJVQf08LRiFrcu3eP27dv8+jRowXbHz16VPVqE3GwfPly1q9fn3XVaRQJgKWOZ8b8qalx4APfmMq6uAwEnbSSyrW3X6b70Ic8MftrvuK1rB0Fx5peY2J+OQOzx+g+9GFJzuPGjRsNd4HIhmnhaEQtVq5cycqVKxdtb7SmVpWk1KmqOyLyVHiDqn4TgT11R3ja6s7c8qy5Hknp62EYhpGPUh3HPwc+EZHf+vkcT+V9RwWoZCOnKAmcx8DssZy5HuXEPbZs2RKZvbWOaeEwLRymRXSU6jiuqmoX8AbwDfDnIvLb6MwqjFrKBM2V65Gtr0cxNa6agyCKYVqEMC0cpkV0lOo43hORF1R1XFWPqeqrqvqTSC0rgFpcdldImZJSalxdvHixIvbWIqaFw7RwmBbRUZLjUNVvVPWTqI1pFPKVKYHq1bgyDMMolqiKHD4uIq9Hsa9iqOWldbmmrkqJe6xdu7ai9tYSpoXDtHCYFtFRsOMIZ4pnoAv4cfnmFEc9VP4sJNO8EOexcePGitpZS5gWDtPCYVpERzEjjqFsT/ixjlcjsKco6qUOTxTO48yZMxW3s1YwLRymhcO0iI5iHMczcZRObxQqueLKMAwjSopxHHtV9Tci8mJceRvpxN3IKWrKWXHV1tZWdXuTimnhMC0cpkV0lFqr6nG8+lQfqeq9yK0qkIGBAa3HfgPhkUR6navwxxVVnSvDMBqHKGpVlRQc95fjfgBsirM6br02ZilmxZX1NF/M2bNn4zYhMZgWDtMiOoqZ63lLRFaISLeIPOU7jI1Av196pOrTV/XemKWQuIfleyymFhNDK4Vp4TAtoqMYx3EQuAtcBT7G6/j3Bl7L2HG8qSsjYvIVSbTAuWEY1abgGIeIvA0cBgZwZdWjM0Skz9/3BNADHFfVsVzv6e/v15GRkSjNSCy54h6wMJAelGiHxox7zMzM1EWOTxSYFg7TwqOqMQ5VPeTHNj72Dx51bGO3qh5V1eOq+g7eCCcnMzMzEZuQXAotkpipRHujjT7Gx62LcYBp4TAtoqPUWlXjqvpJxEtzD4hIRzFvSO/w1QjkW7KbHjgfbtkHNFbs4+bNm3GbkBhMC4dpER3FTFX9TFV/GXrcDXTgTS/1AwfLWZorIvvxstODkcZ7qno3y+v2A6xevbr/3XffBaCnp4cVK1Zw7tw5AFatWsW2bdv49NNPAWhpaWHnzp2Mjo5y755n5sDAALdu3eL69esAbN68mdbWVs6fP4+/f5588klOnz4NeCVOduzYwfDwcCprfXBwkK+//pobN24AXs3/5ubmVCXOtWvXsnHjxlTWaltbG4ODg5w9ezYVrNuxYwfj4+OpE3vr1q3Mzc1x6dIlANatW8f69etTq0La29vZc/xWSpPw1JWnUfZlu//+n3l9l1euXElfXx+nT59mdnYWgGeffZYLFy5w584dAHp7e7l//z5jY96MYXd3N11dXYyOjgLQ2dlJb28vp06dQlUREXbt2sW5c+eYnJwEoK+vj4mJCa5du1a1z+nBgwe8/PLLificBgYGOHPmTGp0vHPnTi5fvszt27cB2L59OzMzM1y5cgWADRs2sGbNGoJl5uV+TlNTU2zYsCGRn1O1v08jIyO0t7cn8nOC6n2fnn/++bKnqopxHPPAJJ6zABC8YPkEMIbXc/ytsowReRMI+pbvzRfjeOqpp/SLL74o55A1T7iveVPTwikriL63eS1w+/ZtVq9eHbcZicC0cJgWHlWNceCtotoNdKlqs6o2qWqXqj6hqj+Owmmo6juq2g8cAT7K955SkhfrjVwl2qExcz7mAhEM0yKEaREdRS3HVdXPK9FbXER2A6PBY1U9Chz3V1pl5eHDh1GbUpNce/vl1BSU5XyQmjoyTIswpkV05HUcfiyDYhxG8J4imAAWOQlVHc3wWiMLlvNhGEY1yOs4VPWaiPy8UGcgIvtwcZCC8B3EmIjs929vAu/me9/SpUuLOUxds27dOsA5j4HZY4umrsptT1srBFoYpkUY0yI6igmO78MbFVzFm1YKAtcdeAl7fww8Tpmrq4qhr69PgxUJjc709PSi6p+BMxhu2UdX04MFo45sq67qIWieSYtGxbRwmBYe1U4APKaq/wqv3Eg/XrmRQ8A/919yWFX/VTWr5dZrkcNSyFTALTz6mJj3YiDhkUe9Tl1ZMTuHaeEwLaKj6ARAP0D+Cz+T/A3/728qETQ3yifsPNLjHmAJg4ZhFE9Nd0JqDme9NTjt7e1Zn8sV94D6K1eSS4tGw7RwmBbRUVIjp6RQr42cKkk4YTA92xysWKJh1DvVTgBMHBbjcAQlGPKRq9YV1Mfoo1AtGgHTwmFaREdNO456b+RUDMVUCs6VbV4PxRIbqWpyPkwLh2kRHSU5DhF5T0ReCT1+Mdxa1kg+mcq01+PowzCM6CkpxiEiLwZ9OULbXlDVTyKzrAAaqZFTPmZnZ2lpaSnpvdmaREFtxj7K0aLeMC0cpoVHnDGOsSCTXET+RkQ+A14qx5BSsKGn4/LlyyW/t95GH+VoUW+YFg7TIjoKdhwickVE/sIfWYwDHf501ZCqPlNuddxSaMRGTtkI+gaUQ73EPqLQol4wLRymRXQUM+L4ADgBvCoiV4B38Bov3amEYUY81NvowzCM6Ck5j0NENuL153gJr4bVR35JkqphjZwcv//97/ne974X6T6LjX3MzyejWVQltKhVTAuHaeERax6H33f8mKq+qqpP4I1AqkotJy9GTSXiPcWOPpIyfWWxL4dp4TAtoiOyPA4/7lFVrJGTI+iFXAmKiX0kYfqqklrUGqaFw7SIjppOADSqR77RRy0Fzw3DKI+adhzWyMmxYcOGqhwnffSRXnE3vVVtHKOPamlRC5gWDtMiOmq6yKE1cnJMTU1VvfpnevC8qWmh04B4Egfj0CKpmBYO08LDihxakcMUcVQJTp++ytWqtpqjD6uY7DAtHKZFdBTtOETkKxE5XGgPcqP+yRf/sNiHYdQXpYw4XgIEOOFnk/8srgKH1sjJsXJl/DUmkzL6SIIWScG0cJgW0VFWjCOUBLjf33QEeK9afcetkVNyCTeMSkrswzCMmGMcIvIUnsM4CEwCb+ONRD4RkX9djlGFMjU1VY3D1ASnT5+O24QFxDn6SJoWcWJaOEyL6CglxvFzEfkKeB+vTlW/qv5YVT/wM8kHgD+P2tBM1PKKsKiZnZ2N24RFxBX7SKIWcWFaOEyL6ChlxLEJ2Kuqm1X1l6r6TYbXvF2mXUYdUe7owwLohpEsio5xZGrY5Mc6Nla7kZPFOBzz8/M0NSV/dXUxsQ9YGP+AwmIgtaJFNTAtHKaFR1wxjiMZtk1k2V5Rpqenq33IxHLhwoW4TSiIYkYfwQhkbOlrC0Yh+agVLaqBaeEwLaKj4D6KIvK6f7crdD/gCWBVOYaIyCTQ4T+8698/qKpZq+7anKXjzp3aaYsSOI/uQx/yxOyv+Qo3+kivuhs4k+Zm+IrXeGL21ynnkW30UUtaVBrTwmFaREcxI44B/9YZuh/cAPaWaoSIdODFTcS/dQIHcjkNo/bJNPoopHCiNY4yjHgpJcbxnqq+GqkRIh2qejf0eA8wqqpjud739NNP6+effx6lKTXL5OQknZ2dcZtRMsXUvQruZ2scVetaRIlp4TAtPKKIceR1HH6+xliQ1JcrSzyKxD9/9PGqqh7N8vx+/ITDtWvX9v/lX/4lAD09PaxYsYJz584BsGrVKrZt28ann34KQEtLCzt37mR0dJR79zwzBwYGuHXrFtevXwdg8+bNtLa2cv78eQBWr17Nk08+mVr/3drayo4dOxgeHk7lkAwODvL1119z48YNALZs2UJzczMXL17Et5GNGzdy5swZANra2hgcHOTs2bOpGM2OHTsYHx/n5s2bAGzdupW5uTkuXboEwLp161i/fj1nz54FoL29nYGBAc6cOZNqTvOHf/iHTE9Pp/oqb9++nZmZmVQPgg0bNrBmzZpUvZ6VK1fS19fH6dOnU1N+zz77LBcuXEgN6Xt7e7l//z5jY57/7u7upquri6CwZGdnJ729vZw6dQpVRUTYtWsX586dY3JyEoC+vj4mJia4du1aQZ/Tn/21qz+Wr+tgcD8cQP/3/2w58/PzvPDCC4n8nHbu3Mnly5er9jl99913rFmzJvLPqRa/T19++SVLly5N5OcElfk+Zfqcnn/++ao4jnlgd7BiSkQmAMVL9gujqlpWnMPf/xBwODwCycaWLVs0+NI2OidPnuS5556L24zICEYgwy376Gp6kHIakNmBzM0lo21t0qi386IcTAuPqqyqUtWm8DJbVe1S1VX+3/CtbKfhs7sQp2HUN8GFf2D2WFHxDyueaBiVp5TM8QkR+T9EpDdqY0RkN97S3oJobW2N2oSapbu7O24TIidT2fYkNo5KMvV4XpSKaREdpQTHe4A/BX4KbASOAu+q6hdlGyPyJvCMqha0QsuC44579+7VffXPYhtHpQfPofGmsBrhvCgU08IjlgRAVR1T1V/4B+4BxoAhEfl9OYb43PX3VxDffvttBIesDxqhE2KxjaOC6atiEwjriUY4LwrFtIiOghMA0/FXV+3F68/xDHCiXGOyraQyjDCFJhCGHUmmBMLwvgzDKJxSYhyvi8gwMI7nNN71g+OR5nYUQktLyX6v7mjE9en5EgiD4Hl6AL2RguiNeF5kw7SIjpISAIEjqvpxZUwqHCtyaASkxz/y5X+ALeE1GpOqxDhE5Km0pL/Xgc9EZGX6rRxDSuH+/fvVPmRiOXXqVNwmxEp6/CPb6qt8S3jrbQTS6OdFGNMiOgqZqhrF1aMCuIY3TXUt7TYepWFGcVhTK4/0/I9MDiTTEt56DaDbeeEwLaIjb5BAVZvSHndVzhyjVETSE/kbl//w37eza9cuug99mCpFEizfhcWB81wVeKG2p7DsvHCYFtFhjZyMuid9BJEvBhL8zVZE0TBqGWvkZI2cUgRFzozFWgTxj1z9zzM1kaqHGIidFw7TIjoS08ipFKyRkyOonmnk1uLa2y+n8j9gYRHFsPMA93hV8wPGml5LVeHN10gqSdh54TAtoiMRjZwMo5qERx/hIoqwOIDeKEF0wyiGRDRyKhWrVeWwOjyOYrVIL+EOi2tgQW2WcbfzwmFaeMQW4/CbO4UNeV1EXijHkFKYC34mGkxMFFxUuO4pVov00UchI5BsbWyTNgqx88JhWkRHKY7j3+EVIwwzCgyVb05xBN26DFIdwYzStMgVQM8VRM9WxiQpDsTOC4dpER2lOI5NqnotvEFVR4G+SCwyjJhJr4GVKZEQFq/CyhQDSYoDMYwoKaVK4AkR+R9V9T8GG0TkFaDqwQZr5OTo6emJ24TEEIUW4ZhFpkTCXKuwmpthrOm1VB5InKuw7LxwmBbRUWojp2HgM7zeGZuAfqA/fSRSaSw47picnLTqnz6V0iJXIUXIXExRldQy3oBqOhA7LxymhUecjZy68Dr/jeFVyl1VbacB1sgpjCU3OSqlRaExEMi/jLda01h2XjhMi+gouaGFqn4QpSF34vk3AAAZHElEQVSGUSukJxFCcbWwmpqoyWRCwwgo2nGIyASQcX5LVauaPW6NnByrVlU9cT+xVEOL9BhIuhNpbl48+lAlNb1VrWx0Oy8cpkV0lBLj2Ji+CTgA/DdV/WVUhhWCFTl0zM/P09RUyiK5+iNOLQIHEA6iB4RHIulxkEoVVLTzwmFaeMQV4xhPu42p6kHgp+UYUgrWyMnx6aefxm1CYohTi/SlvJmW8cLicibZCiqWGwex88JhWkRHKVNV3Rk2d+KtrDKMhidwHrmW8Qakx0FWNT9IxUHCy3ktBmIkiVKCBKN4MY5wV5Q7wMFILCoCa8zisHiPIylahB1IMA2VKYge3IeFcZAoGkslRYskYFpER9ExjiRhMQ6jlkifdgqXdA/IFAcJ/saZD2LUD1HEOGracWzdulUvXrwYtxmJYHR0lL4+q/oCtaFFejJhpkA6LJzWCh6nB9IhuxOpBS2qhWnhEYXjyDt2y7X8Np1qL8e16riOe/fuxW1CYqgFLXJNY2WawoL8+SCZnEctaFEtTIvoKGTSL7z8VoBDwFXgRGjbAeC/RWuaYdQ/mRxIelfCgGz5IOXGQQyjWErJ47iiqpsL3V7kvveEH6vq8Vyv7+vr09HR0XIOWTdMTU3R3t4etxmJoJa1yDSFBfnzQYK/6dNY5//XXTWrRdTU8nkRJVWZqsrAqgzVcV+gzJ7jIvImMKaqx0WkA/gYyOk4Hj16VM4h64pbt27Zl8KnlrXINAKB3KOQcHmToDJvMIW1/X8/tWjfjUotnxdJoxTH8SrwNyLyEa467ouU33P8LVXtBFDVuxSQF/Ldd9+Vecj64fr162zatCluMxJBPWhRaFl3yBwHCUqaAIko754E6uG8SAolraryRwR7gQ68boDvqeo3JRshshuvg+Bhf399wHFVHcvw2v3AfoDVq1f3v/vuu4BXa3/FihWpCpirVq1i27ZtqWzRlpYWdu7cyejoaCpINjAwwK1bt7h+/ToAmzdvprW1lfPnz+PvnyeffJLTp08DXv+PHTt2MDw8zNTUFACDg4N8/fXX3LhxA4AtW7bQ3NxMsNpr7dq1bNy4kTNnzgDQ1tbG4OAgZ8+eZXp6GoAdO3YwPj7OzZs3Adi6dStzc3NcunQJgHXr1rF+/XrOnj0LQHt7OwMDA5w5cybVBVFVWbNmDbdv3wZg+/btzMzMcOXKFQA2bNjAmjVrCJYvr1y5kr6+Pk6fPs3s7CwAzz77LBcuXODOnTsA9Pb2cv/+fcbGvI+hu7ubrq4ugunBzs5Oent7OXXqFKqKiLBr1y7OnTvH5OQkAH19fUxMTKS6r1Xjc3rw4AEvv/xyIj+nnTt3cvny5ZI+p//5/3VfsUwrsaC41Vj/9/+0PtbPqdrfp5GREdrb2yv+OSX9+/T888/Xx3Jc3xkcUVXxH3cAI6qa8+dBb2+vWqlkjxs3brBu3bq4zUgE9a5FOA4STGFB5uW8kL8uFjTGKKTez4tCiSXG4cczDrFwKkkALWM57ph/A29Hd0WkR0R6Mo06QraUeLj6w7ohOupdi0wlTSD3NFY4FhLEQaCxprHq/byoJqWsqrqD5zhOpD+nquMlGeF1FfwoPMIQkUm8roJZHceWLVs0mCZodE6ePMlzzz0XtxmJoBG1yDQKyfS7Kj24Dpm7FEL9OZFGPC8yEdeqqklVPZb/ZYWjqmMicjd47E9VjeVyGoZhODKNQgIHMj9P3j7p4WC6NZky8lGK4zgoIn8BvAcsGGGU2T52r4gM4SUXbqKAVVpLliwp43D1xerVq+M2ITE0shbZprFgcYMpWFydF5wTSZ/GCu+/Fmnk8yJqSpmqmsjyVDkxjpLo7+/XkZGRah4ysczOzlr1Tx/TwhFoEVz8SwmmQ31MZdl54RFXI6euLLeq92UMlvAZpJY4GqZFmECL4CI/MHuMnu8yN5kKbulNplJTWUtfW9Rsqpaw8yI6Sna/IvIU0ANcVVVbE2sYCSdfUiHkTyyE7PWx0o9h1C+lLMd9HK8cyON4MY4efwXUi6pa1fKT1j/YYUsNHaaFI5sWpZQ2CQiKLAZxkPQqveH9Jwk7L6KjlBjHu8CJ8Moqv87Ui6r6k4jty4k1cjKM6CikwGLwONNlI1NyYRIdSKMTSyMnEbmTKZ6RbXslsUZOjuHhYQYGyjoX6gbTwlGKFunTT8VU6YXs2ekQryOx88IjrjyOcRF5XlX/NmTIC6Qtza0G1sjJYQsFHKaFoxQt0mMhYQeQKTs9IFN2ejCVlYQMdTsvoqPU6rjDGarjvhilYYZhxE82J5LJgYSD6mEHAgsdiQXVa5+SixyKyD48p3GVMqvjloo1cnJMT0/T1tYWtxmJwLRwVEKLTLGQcHZ6QL6prPBIJKCSDsTOC4+4YhwvABOq+kVo2+t4JUI+KceYYvnhD3+oX375ZTUPmViuXLnC5s1lNWCsG0wLRyW1yJTHka3cO8QfVLfzwiOuGMcR4KW0baP+9mfKMaZYrJGT48aNG/al8DEtHJXUIv3inmkqK71OFiweiaTHRCo1lWXnRXSU4jg2pdekUtVREemLxiTDMGqRbLkhUNtBdWMxpUxV/Q3wF2k9x18B/rzc4U+xWCMnxz/+4z/yB3/wB3GbkQhMC0cStMgUEyl3KiugGEeSBC2SQFxTVW/grap6A7eqqp8CeoRHjTVycjQHy1cM0yJEErTINBLJVvId8k9lZcpUDx8nG0nQol4o2nH4PTK6RORP8WpVnVDVDyK3rACCHsMGXLx40cpG+5gWjiRpka/ke6FTWVBa/5AkaVHrlFzkMC5nYRhGbVNIUD1ffkih/UMyHc8on5ouTm+NnBxr166N24TEYFo4akGLYqeyAtKbUIV7qcPiwPpf/9kmjGgoOQEwCVgjJ8fMzIxV//QxLRy1qEUp+SGQPbBeD02ooiSWBMAksWXLFr106VLcZiSCkydP8txzz8VtRiIwLRz1oEUxmephMnUyzLQ6q9GcSFVWVYnIykJ3Vu1+HIZh1D+58kOy9Q+BzP3Um5vhK17LmCOSfjwjO4XEOK4BCmRb+xo8p0BVy6pbIyeH1eBxmBaOetIiU1A9UxdDyN3JMFPhxVKX+DYqNT1VZY2cDMMoJiaSrfBiQLCtnuMiDR/j+MEPfqC/+93v4jYjEZw9e5bBwcG4zUgEpoWjEbUIO5JgKgtyZ6tnciT1GheJK3McEXkKGAA6wttV9ZflGFMs8/Pz1TxcorFkSIdp4WhELbIlGpayxBcWLvO1uIhH0Y5DRH4OvAWcwGve9D6wG68vR1Udh2EYRjbSm1AVkq2eLS6S7kQaPS5SUs9xYKOq3hORYVUdEJEO4F1V/UlFrMyC5XE4anG9fqUwLRymhSOsRaFLfGGxE0mn1uIicTVymlDVLv/+e8C/VdVzInJHVau6qmr79u16/vz5ah4ysfzDP/wDf/RHfxS3GYnAtHCYFo5MWmQKrNd7XCSuGMdREXlFVX8DvA38rT8KGS/HEBEZAt4E7gLDwAG/oGJWHj16VM4h64qbN2/aBcLHtHCYFo5MWuRa4gsLnQhYXCSglOq4h0L3R0VkAG/q6uMybflMVa1OumEYsWFxkcJIzHJcEdmjqseLec9TTz2lX3zxRf4XNgC3b9+2ktE+poXDtHCUq0W14iJQWUcSy1SViHQDB4BFrWLLDI73iMgevKmql4DDqno31xuS4vSSwNzcXNwmJAbTwmFaOMrVopDSJ5C/p0i2kvDZyqAkcTRSSnD8CvAx8FH6c+X06BCRjsBR+P3Lh1T1pQyv2w/sB1i9enX/u+++C0BPTw8rVqwgaCW7atUqtm3bxqeffgpAS0sLO3fuZHR0lHv3vJJaAwMD3Lp1i+vXrwOwefNmWltbCQLuq1ev5sknn+T06dMAtLa2smPHDoaHh5mamgJgcHCQr7/+mhs3bgCwZcsWmpubuXjxIuCVtd64cSNnzpwBvBIQg4ODnD17NrXGfseOHYyPj3Pz5k0Atm7dytzcHEEBx3Xr1rF+/XrOnj0LQHt7OwMDA5w5c4aZmZlAe9asWcPt27cB2L59OzMzM1y5cgWADRs2sGbNGoJM+5UrV9LX18fp06eZnZ0F4Nlnn+XChQvcuXMHgN7eXu7fv8/YmBdq6u7upquri9HRUQA6Ozvp7e3l1KlTqCoiwq5duzh37hyTk5MA9PX1MTExwbVr16r2OT148ICXX345kZ/Tzp07uXz5ctU+p6mpKTZs2JDIz6na36eRkRHa29sj/5z+7K9dDAQyx0XCFDoayRRg/y9//qOyP6fnn38+3lVVlURENF/Mw6rjOuqhCmpUmBYO08JRDS1ylT+Zn3d1sgLSA+vZVmlFudQ3rlVVh0XkfwHej6oarj/COKaqRfUtX7p0aRSHrwvWrVsXtwmJwbRwmBaOamiRq7shuMq8kL1eVqZVWquaH+St6pvp+JWilBHHn+JliwdvvItfHbfUPA4/gfBVVT3qP96Ntxx3b6739fX1aTBt0uhMT0/XVSXUcjAtHKaFI24t0i/0heSLBGQbiQQUE2SPKwHwDnAIr+TIRPg5Vf2mZEM8Z9HjP9xEAcFxm6py2JSEw7RwmBaOJGmRzYnkWqWVq6oveNvTg+xhAicS11TVpKoudm1loqonot6nYRhGEik0XwTy54wE24O+65A7AXHp2ieKCglkohTHMSQifwG8R1q2uKpeK9egYmhOjzQ1MO3t7XGbkBhMC4dp4UiqFnnjIhlyRrLFRjIt+U1PQOyj/IrJJa2qyvJUyTGOUrFGToZh1DP5amkFFLPk95ljUwz/17myqnSUUnKk4ktxC+XBgwf5X9QgnDlzhh07dsRtRiIwLRymhaMWtchXSysgfclvrgTEKCglc/x1Vf1VNIcvD2vk5AgSlwzTIoxp4agHLTKtkkqf2home2HGqIptlBLjeENErqrq30ZjgmEYhlEqhQbaMyUglkopMY6n8QLjb+OVHknFPKJKCCwUa+TkmJ2dpaWlpE7AdYdp4TAtHI2oRXqM5KuW1/jv/q/yYxxNJbznY2AV8AtgFG9l1TXK7MdRCvUw9IyKy5cvx21CYjAtHKaFoxG1uPb2y6kbeMtyv9SePO/KT9GOQ1W70m6rgr9lW1Mk1sjJERRjM0yLMKaFo9G1CBzIdze/KnuappQRByLyuoj8n6HHr4jIU+UaYxiGYSSfoh2HiBwGnsDrmREwDkSeTZ4Pq8Hj2L59e9wmJAbTwmFaOEyL6CglUrRfVVf5taUAUNXPRaT8ibMisUZODov3OEwLh2nhMC2io5SpqnG/C2Dqqi0irxBDcPzhw4fVPmRiCRrMGKZFGNPCYVpERykjjn3AcWCjiPwM+DHQ798MwzCMOqeUVVWf+yV53wK+BxzxV1Zdi9q4fFgjJ8eGDRviNiExmBYO08JhWkRHydkwQdOlOFmyZEncJiSGNWvWxG1CYjAtHKaFw7SIjlJWVb0gIn8jIndCtwm/wVNVsSKHDqsS7DAtHKaFw7SIjlJGHO/jdQA8ELEthmEYRg1QiuMYr0QHwFKwRk6OlStXxm1CYjAtHKaFw7SIjrxFDkUkXe2fAn0koAOgNXIyDMMojih6jhcS47iGK2R4DRjCcx4f4BU5DG5VL1M7NTVV7UMmltOnT8dtQmIwLRymhcO0iI68U1VJ6viXjmWOO2ZnZ+M2ITGYFg7TwmFaREfeEYeIfFYNQwzDMIzaoJAYx0RSRx0W43DMz8/T1FRSseO6w7RwmBYO08KjWjGOxM4HTU9Px21CYrhw4ULcJiQG08JhWjhMi+goZDlup4g8D+RtNaiqn5RvUuHYnKXjzp2q518mFtPCYVo4TIvoKDSP44MCXqN4LWUNwzCMOqYQx3G32jEOETmiqnkz05ctW1YNc2qC3t7euE1IDKaFw7RwmBbRkbgYh4j0Aa8W8tq5ubkKW1M73L9/P24TEoNp4TAtHKZFdBTiOPLGNqJCRDqKeb119HKMjY3FbUJiMC0cpoXDtIiOQhzH4Ypb4ditqqNVPJ5hGIZRJHnzOKqF38N8WFXvisikqnZmed1+YL//cDtwvlo2JpzvAb+P24iEYFo4TAuHaeGxRVVXlLODRDgOf4qqJxht5HIcae8bLjeRpV4wLRymhcO0cJgWHlHoUHIHwIjZDSAiPf7jDn9kcUJVbWLSMAwjQSTCcajq8fBjEUlEa1rDMAxjMYkq3CIiHSLypn//zdAIJBvmXBymhcO0cJgWDtPCo2wdEhHjMAzDMGqHRI04DMMwjORjjsMwDMMoCnMcNYSIHEl73OPHgnb7fzsKec4wDKMcEhnj8OtVdQEdwEvAULAs1w+Y78Hrc94HHFXVu/meq3V8TT4O57f4DVn6/fsdwDFV3ZvvuVpFRIaAN4G7wDBwoJHPCxHZE34crE5sNC1EZBLvWgHeudEBHFTVdxpQiz5gAJgAeoDjFfmOqGribsAk0OHf3w98FHpuJHS/A3i/kOdq+eb/L33AZGhbX1iXQLd8z9XyDdiT47mGOi/wHOie0P8U/h8bRgv/f9idtm1/I2oRnBdpj49UQoukTlVtVOftJvBGH4E3nQhe5L9md77n6oBMNbx68H5dhZnwfznkeq7uaNDz4i31Rxiqelfd6LLhtFDVE8F9fxR2wr/fcFoABzJNS0etRSIdhy4cIh0ADvr3G+5i6dfwOpHhqWw9UjryPFfL9IjIHj9uMxT6gjTUeeGfE2MhLcI5Tw2lRfha4Z8PXeqqTTSUFj5DwLiI7Perb1Tk2pmIzPFMhObcPgr9omioi6X/RZjQzHONEyz+37oKeK6WCc/JTgDv48XAGuq8wPui96mLaQwDI8AmGk+LMG+xsJp3w2mhqkf960bQCO8EnlOIVItEjjgAVHVMVd8B7orIR/7mRrtY7sb9yt6DX8PLd6oZa3j5U1q5nqtZwg7U/1+C4XSjnRdjhD5jX5ce/7xoNC3C7E77kdVwWojIm6r6jj91eQSoyLUzcSOOYKThOw2A94Aj+S6WIpn7TdXyxVLz1PASka7Q/R78KS1fj4zP1Sr+POyxYC4/jYY6L8j8/97N8Vw9awGkpu8m0jY3lBa+Bin7/dHHJv+7E6kWiXMceMPwVWmP76pbUtYwF8sAf+i537//Jm6J3T7/cbCEbl/obbmeq0XG8H5BAakvSRAczvrZ1+N5oapjIpI+tz/WyN8RvHN8wZRuo50XeI5zUUxUXbuKyLRIah7HHtxQ6SXgcOif78MTJ9Na5KzPGbWP7yyCgN0mvPMi72dfj+eF/+U+AFzF0+JIyHE0lBaQavC2SVUPpm1vKC3Srp0deK0pIr92JtJxGIZhGMklscFxwzAMI5mY4zAMwzCKwhyHYRiGURTmOAzDMIyiMMdhGIZhFIU5DsMwDKMozHEYseAX5psM3dS/BY9H4raxHimniF+mqqtGY2KOw4gFv55OZ3DDy/p9KbQtU2kRowz8OkYZS08UyKvmPAwwx2EYZSMiR/zM5Woc6yO/pPx+/6/6998Ukfclrb1w6H378Uu0pG/393HEzx4Otu/2R35DwTa/Ttpblfi/jNoiibWqDMPIgH9hHwraDATOIK3wZTYH1h9+XYj38MqVBGW4gympDg21KQ5xVUR6yhy5GDWOjTiMxCIifSIykv7L13+8X0Su+vf3hGImV8Pz+P62N8OvLXD/wXNBccmhUPzlzdBr38crQDnkH2N3aB89/v0e8fpi59p/RlvS6Al3uwP24spmBwxn0hGvplUmFhTF81/bk16ZOcQJvD45RgNjjsNIMseAvcEv37Rf0/2qugmvw9n7eBWUO/EubAfSd+S/dh/wfmiePtf+3wcOhn6lf+a/biPwVjCto6p7gaP+azelXdhzkb7/XLYE/0P6xTxrJdQMr8tWIvslfOfjO6y7ucpp+yONl7I9bzQGNlVlJBL/wtmDd6EHv9Jn6CXv+3+Dbe/5fz9iseMIyq8fF5FRYLdfRjrX/ofCTiDc31tETgADZL8YF0L6lFMuWxbhO66xAqeMVpGlHwPwKl4J/t14I4nPcrzWMABzHEZyCcqmv5Pl+eDiNgGL+tTnYhjX8yXX/hdcuP0L9Vv+e3vwLrDlEN5/vv81E9l60Wci40oofyot6OUx6o843iJDEN0wwthUlZFUPqPwKZFieigM4DmdfPtPdZPzL7Af413c+3Gjm3IId6sr5n8NSE0xFcBdMjuP3XjB9VFIrZrqCa+uMoxMmOMwEok/NTQQBLP94HGpQdlgH/vxA79F7r8HmPB/lXfgep0H3AU2iUhHKDAf7uO8KOYSpsT/tZgRx2e4BlhhMgXXD5N/yW1NNzsyyscch5Fk+vEC0ZN4weNS595XichVvED6i8Xu349F3PVf936G132Et7JqBHeBPgKM+BnwdwqwsSBb/NVkQXxnf4GjgxPAM6F9dPgrw3YD/WlJfZuAPeGVY2nH3w28W8AxjTrGOgAadY1/Ie5v9LwDEVmQqxH3fozaxkYchtEYHCw3u92fhis0rmLUMeY4DKMB8FedDZdT5BDoy5EYaDQQNlVlGIZhFIWNOAzDMIyiMMdhGIZhFIU5DsMwDKMozHEYhmEYRWGOwzAMwygKcxyGYRhGUfz/TF0JxyCs5LIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "j = problem.fluxes[\"HeatFlux\"].function\n",
    "g = problem.gradients[\"TemperatureGradient\"].function\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
397
    "k_gauss = -j.vector().get_local()[::2]/g.vector().get_local()[::2]\n",
398
    "T_gauss = problem.state_variables[\"external\"][\"Temperature\"].function.vector().get_local()\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
399 400
    "A = material.get_parameter(\"A\");\n",
    "B = material.get_parameter(\"B\");\n",
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
    "k_ref = 1/(A + B*T_gauss)\n",
    "plt.plot(T_gauss, k_gauss, 'o', label=\"FE\")\n",
    "plt.plot(T_gauss, k_ref, '.', label=\"ref\")\n",
    "plt.xlabel(r\"Temperature $T\\: (K)$\")\n",
    "plt.ylabel(r\"Thermal conductivity $k\\: (W.m^{-1}.K^{-1})$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "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",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
434
   "version": "3.6.9"
435 436 437 438 439
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}