tutorials/tutorial05/05_SPDE_on_fenics_solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 05 SPDE on fenics\n",
    "\n",
    "We consider the following stochastic Poisson problem on the square $\\mathbf{x}=(x, y)\\in\\Omega := [0,\n",
    "1]^2$:\n",
    "\n",
    "$$\n",
    "\\label{eq:sPDE}\n",
    "\\begin{cases}\n",
    "-\\nabla\\cdot (\\kappa\\ \\nabla u)=1, & \\mathbf{x} \\in\\Omega, \\\\\n",
    "u = 0, & \\mathbf{x}\\in\\partial\\Omega_{\\text{top}}\\cup\\partial\\Omega_{\\text{bottom}},\\\\\n",
    "u = 10\\ y\\ (1-y), &\\mathbf{x} \\in\\partial\\Omega_{\\text{left}}=\\partial\\Omega_{\\text{inf}},\\\\\n",
    "    \\mathbf{n}\\cdot\\nabla u = 0, & \\mathbf{x} \\in\\partial\\Omega_{\\text{right}}=\\partial\\Omega_{\\text{out}},\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "with homegeneous Neumann boundary condition on the right side of the domain (outflow), that is $\\partial\\Omega_{\\text{right}}$, Neumann boundary conditions on the left side of the domain (inflow), that is $\\partial\\Omega_{\\text{left}}$, and Dirichlet boundary conditions on the remaining part of $\\partial\\Omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diffusion coefficient $\\kappa:(\\Omega, \\mathcal{A}, P)\\times\\Omega\\rightarrow \\mathbb{R}$, where\n",
    "$\\mathcal{A}$ is a $\\sigma$-algebra, is such that $\\log(\\kappa)$ is a Gaussian random field, with covariance function $C(\\mathbf{x},\\mathbf{y})$ defined by \n",
    "$$\n",
    "C(\\mathbf{x}, \\mathbf{y}) = \\exp\\left(-\\frac{\\lVert \\mathbf{x} - \\mathbf{y} \\rVert^{2}}{\\beta^{2}}\n",
    "\\right),\\quad \\forall \\mathbf{x},\\mathbf{y}\\in\\Omega, \n",
    "$$\n",
    "where $\\beta=0.03$ is the correlation length. This random field is approximated with the truncated Karhunen–Loève decomposition\n",
    "$$\n",
    "\\kappa(s, \\mathbf{x}) \\approx \\exp\\left(\\sum_{i=0}^m X_i(s)\n",
    "  \\gamma_{i} \\boldsymbol{\\psi}_i (\\mathbf{x}) \\right)\\qquad\\forall (s,\n",
    "\\mathbf{x}) \\in \\Omega\\times\\Omega, \n",
    "$$\n",
    " where $(X_{i})_{i\\in 1,\\dots, m}$ are independent standard normal\n",
    " distributed random variables, and $(\\gamma_{i},\n",
    " \\boldsymbol{\\psi}_{i})_{i\\in 1,\\dots, d}$ are the eigenpairs of the\n",
    " Karhunen–Loève decomposition of the zero-mean random field $\\kappa$.\n",
    "\n",
    "In our simulation the domain $\\Omega$ is discretized with a triangular\n",
    "unstructured mesh $\\mathcal{T}$ with $3194$ triangles. The parameter\n",
    "space has dimension $m=10$. The simulations are carried out with the finite\n",
    "element method (FEM) with polynomial order one, and for each\n",
    "simulation the parameters $(X_i)_{i=1,\\dots m}$ are sampled from a\n",
    "standard normal distribution. \n",
    "\n",
    "The solution $u$ is evaluated at\n",
    "$d=1668$ degrees of freedom, thus $(V, R_{V})\\approx(\\mathbb{R}^{d},\n",
    "S+M)$ where the metric $R_{V}$ is approximated with the sum of the\n",
    "stiffness matrix $S\\in\\mathbb{R}^{d}\\times\\mathbb{R}^{d}$ and the\n",
    "mass matrix $M\\in\\mathbb{R}^{d}\\times\\mathbb{R}^{d}$. This sum is a\n",
    "discretization of the norm of the Sobolev space $H^{1}(\\Omega)$. The number\n",
    "of different independent simulations is $M=1000$. \n",
    "\n",
    "This test case was remodelled from [Zahm, Olivier, Paul G. Constantine, Clementine Prieur, and Youssef M. Marzouk. \"Gradient-based dimension reduction of multivariate vector-valued functions.\"](https://epubs.siam.org/doi/pdf/10.1137/18M1221837).\n",
    "\n",
    "The library used for solving the previous PDE is [fenics](https://fenicsproject.org/). Mind that running these simulations could take a lot of time, around 15 minutes on a laptop with a 6 cpu cores processor and 16 GB of RAM. All the additional .npy and .png files will be added in the subdirectory data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sys\n",
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_mesh_map(mesh, dim):\n",
    "  m_map = np.zeros((dim, 2))\n",
    "  for j, cell in enumerate(cells(mesh)):\n",
    "    m_map[j, :] = cell.midpoint().array()[:2]\n",
    "  # print(m_map.shape)\n",
    "  return m_map\n",
    "\n",
    "def compute_cov(mesh, beta, dim, mesh_map):\n",
    "  print(\"start covariance assemble\")\n",
    "  cov = np.zeros((dim, dim))\n",
    "  for i in range(dim):\n",
    "    for j in range(i, dim):\n",
    "      cov[j, i] = cov[i, j] = np.exp(-(np.linalg.norm(mesh_map[i, :]-mesh_map[j, :], 1))/(beta))\n",
    "      #cov[j, i] = cov[i, j] = np.exp(-(np.linalg.norm(mesh_map[i, :]-mesh_map[j, :]))**2/(beta)**2)\n",
    "  print(\"end covariance assemble\")\n",
    "  evals, evecs = np.linalg.eig(cov)\n",
    "  E = (evals[:m]*evecs[:, :m]).T\n",
    "  return cov, E\n",
    "\n",
    "def set_conductivity(sim_index, mesh, c):\n",
    "  # print(\"set conductivity\")\n",
    "  D = FunctionSpace(mesh, \"DG\", 0)\n",
    "  kappa = Function(D)\n",
    "  dm = D.dofmap()\n",
    "  for i, cell in enumerate(cells(mesh)):\n",
    "    kappa.vector()[dm.cell_dofs(cell.index())] = np.exp(c[sim_index, i])\n",
    "  return kappa\n",
    "\n",
    "def boundary(x):\n",
    "  return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS\n",
    "\n",
    "def boundary0(x):\n",
    "  return  x[0] < DOLFIN_EPS  \n",
    "\n",
    "def compute_solution(sim_index, mesh, kappa, pl=False):\n",
    "  # print(\"compute solution\")\n",
    "  V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "  u0 = Expression(\"10*x[1]*(1-x[1])\", degree=0)\n",
    "  bc = DirichletBC(V, Constant(0.0), boundary)\n",
    "  bc0 = DirichletBC(V, u0, boundary0)\n",
    "  \n",
    "  u = TrialFunction(V)\n",
    "  v = TestFunction(V)\n",
    "  f = Constant(1.0)#Expression(\"exp( - 2*pow(x[0]-0.5, 2) - 2*pow(x[1]-0.5, 2) )\", element=V.ufl_element())\n",
    "  a = kappa *inner(grad(u), grad(v))*dx\n",
    "  L = f*v*dx\n",
    "  u = Function(V)\n",
    "  solve(a == L, u,[bc, bc0])\n",
    "\n",
    "  if pl:\n",
    "    u_pl = plot(u, title='u')\n",
    "    plt.colorbar(u_pl)\n",
    "    plt.show()\n",
    "  return u \n",
    "\n",
    "def restrict(mesh, v):\n",
    "  # print(\"restrict on outflow right side\")\n",
    "  Right = AutoSubDomain(lambda x, on_bnd: near(x[0], 1) and on_bnd)\n",
    "  V = FunctionSpace(mesh, 'CG', 1)\n",
    "  bc0 = DirichletBC(V, 1, Right)\n",
    "  u = Function(V)\n",
    "  bc0.apply(u.vector())\n",
    "  v_restriction = v.vector()[u.vector() == 1]\n",
    "  return v_restriction.mean()\n",
    "\n",
    "def compute_gradients(component_index, mesh, kappa, E, boundary, cache, solution, pl=False):\n",
    "  # print(\"compute gradient\")\n",
    "  V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "  bc = DirichletBC(V, Constant(0.0), boundary)\n",
    "  w = TrialFunction(V)\n",
    "  v = TestFunction(V)\n",
    "  a = kappa*inner(grad(w), grad(v))*dx\n",
    "  D = FunctionSpace(mesh, \"DG\", 0)\n",
    "\n",
    "  dkappa = Function(D)\n",
    "  dm = D.dofmap()\n",
    "  for i, cell in enumerate(cells(mesh)):\n",
    "    dkappa.vector()[dm.cell_dofs(cell.index())] = kappa.vector()[dm.cell_dofs(cell.index())]*E[component_index, i]\n",
    "  rhs = dkappa * inner(grad(solution), grad(v))*dx\n",
    "\n",
    "  w = Function(V)\n",
    "  solve(a == rhs, w, bc)\n",
    "\n",
    "  if pl:\n",
    "    w_pl = plot(w, title='w')\n",
    "    plt.colorbar(w_pl)\n",
    "    plt.show()\n",
    "  return w\n",
    "\n",
    "def show_mode(mode, mesh):\n",
    "  c = MeshFunction(\"double\", mesh, 2)\n",
    "  # value = mode.dot(E)\n",
    "\n",
    "  # Iterate over mesh and set values\n",
    "  for i, cell in enumerate(cells(mesh)):\n",
    "    c[cell] = mode[i]#np.exp(value[i])\n",
    "  \n",
    "  plot(c)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start covariance assemble\n",
      "end covariance assemble\n",
      "Karhunen-Loève mode shape (10, 3194)\n"
     ]
    }
   ],
   "source": [
    "# Read mesh from file and create function space\n",
    "mesh = Mesh(\"data/mesh.xml\")\n",
    "\n",
    "#dim = 6668 #mesh_2\n",
    "dim = 3194\n",
    "m = 10\n",
    "M = 500\n",
    "d = 1668\n",
    "cache = np.zeros((d, m))\n",
    "cache_res = np.zeros(m)\n",
    "\n",
    "#choose lengthscale\n",
    "beta = 0.015 #beta=0.03\n",
    "inputs = np.random.multivariate_normal(np.zeros(m), np.eye(m), M)\n",
    "\n",
    "#samples \n",
    "np.save(\"data/inputs\", inputs)\n",
    "\n",
    "#covariance modes assemble\n",
    "m_map = compute_mesh_map(mesh, dim)\n",
    "cov, E = compute_cov(mesh, beta, dim, m_map)\n",
    "c = inputs.dot(E)\n",
    "np.save(\"data/covariance\", cov)\n",
    "np.save(\"data/cov_modes\", E)\n",
    "\n",
    "print(\"Karhunen-Loève mode shape\", E.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mode number 2 of Karhunen-Loève decomposition\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD4CAYAAADMz1tMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABcCElEQVR4nO29Waw0yXUe+J1carvrv/febEokTc1YEuSmJBiaGcmCxqSMAWFAD7KEEUawQRAQDT/4QcI82AP4ZQxjAMPQ0iAEQjAwMDGABQ89blsYzMCWYZkyW2OJIltmq9VUd/+9/ft/l9pyOfOQGVkRkScis+rWXeq/+QGJeyszMjIqK744J06cc4KYGR06dLhcCM67AR06dDh7dMTv0OESoiN+hw6XEB3xO3S4hOiI36HDJUR0Xg++fv06f+xjHzuvx3focCnwB3/wB/eY+YZ9/tyI/7GPfQyvvfbaeT2+Q4dLASJ6WzrfqfodOlxCdMTv0OESoiN+hw6XEB3xO3S4hOiI36HDJUQj8YnoK0R0h4i+5bhORPRPiOhNIvomEf3Q+pvZoUOHdaKNxP8tAJ/1XP8cgE+UxxcA/MbJm9WhQ4fTRCPxmfl3ATzwFPk8gH/KBb4OYJ+Inl5XAzt06LB+rMOB51kA72qfb5fnPrALEtEXUGgFeOGFF1pV/tKv/m8nb+EFANOSN5CQJ0HVsUwKhaUf3NyMDueP737p757o/nUY96SeJXYXZv4yM7/MzC/fuFHzInyisTSBfHy1r52M20uB6cRjSYcLgHUQ/zaA57XPzwF4fw31PlFYiixqlFjqnoa6VoSr3R35NxvrIP7XAPxCad3/UQCPmbmm5l9mrET66rOjXNhwfY1Q7be/R0f+zUXjHJ+I/hmAHwdwnYhuA/j7AGIAYOZXALwK4KcBvAlgDOAXT6uxm4gTkb46D3mIbls38YlZykHZNqseps4OsIloJD4z/42G6wzgl9bWoicEXp7pTFEFw/JcbpX16WSBVZ7gN/wFwkX7lN5w0cBYHwDUv/qlbjC42Di3sNwnDU6i29LWZkTApuRWRJfqkyS3Kq8GgBAOcH1QUc/Rm9SWsUJb9I8d+S82OpfdE8Jr5a6MdLyQtiQctfs8D3RJ4ZDBsUdCA8Wv7fvFXe1Zpi0auhWAi4uO+EtCN3S17tSqnKRqrwJi89DbFzaQH5AHgI6glwqdqr8EXNZtL5Z23GnRDo20lNfPLS42DDTKRrCMz4BdZQvDoXS5mwacLzrit8DK6mrTfZIm0ESitjpaE7PU9RBAbq/Twd32JgPiGprW4fTREb8Ba3G1ra5ZnyXVn7heLidwlIuNkQYCDrnSBFq1LWCT/G0GrCaDYIsX1xkAzw8d8U+KZXquLkmVIa3F7RyVLLYGCorkmynIwYlp3mdmUNpyFLNjAqR2NsUNOKYAa3Ap6LAGdMa9BqwskdrOmxss6ewwCLpIX12Ps3pdvntE7QPmQLUMHOzuSH8x0BG/Aa0dcYzzvntano/zOum5IHwT6asqG8jPweJoV2HLch0uPDpVf1X4wmbF8v7q9GU4KglP1tIcxTmIGHkqe+lQqEz8+j0Z8rml9kdcM+hxANkuUHsITmzc63D+6IjvQaNjjnHOU5Gh2gseb9Lau11FvGBlEC0kedQv/k+m8k8ZxDmCOEd6HBf1qCVJu35G4TacUzOx9eafcBCw3X07nA064q8C2xy9rKQvyd+G8IBJemYClc+OeosBIB6kABYDQBCb4jsoy9pGP9Ucsc1tmmfHC9iVtpzUdxb+s0VHfAca+6sqsKI3nov0ZFvuY1n/1kmvQw0AWSpP3NW8n7Pip/eSrc0AwA3lljDjd+Q/O1xq4jcFojViVdJXgTVUr6OFFuAifXU9zJzE19H6O7vm9R1JNxaXyqqv/OslP/uVlplW6PguCzqF7JTubdsQhRmiMAMzoddLq/N5Yj2UAIqEZ9kefNY93ilNNwhsFJ4oie9TFdsQu33QzYqSXliGsy33NgLheqBpCWGYI0lDRGFdC+j3EwBAGoZIZsJPrQYaNTC0Uct16a+KNo1XLdX9Ts0/OzxxEt+W5msLDXVEw61V0lkSNwhzsEWqMM6QW18ojkzSE3FlAJTAtgYQ58V7CuCN/Ovw5OCJI75CK8Iv07lXGT2sWygTfO1tEpYIwroYDQWHnKpuWizVSYj7qfuink2nFq7b4v20eTXdIHKh8MQSfymclPxtblfr1Rb5Sa2dK2nPQn1qgOLF9TwnhNbgwFwnf6o5+8T9dDHQuIwdOEXyN6Bz5z07dMRXOA3J7/LoFSQ/AEPVz1tY5ZMkRJKEYDZJ75P8CLi2ksCu9uio1uV9ZRznl5hvdVl7zgYd8XWsUx1l6+8KyNPAMP5lmfxzpWloSPY0C5DlxaFjPo4XH8J8cQhwJvZQ1n390MvWEncsz+JuVnD6eGKI36p/nZXByn6E5OGbkVvyqzIhI0+D6tCRC0tvaRYgdQwOWRoi7C23XFitQhhRekuo/UuSvrMnnh2eiOU8o3+1cRU9zaDwZTpuwEBKgCKYPWfPyJD4agDwGfkkZI6gngLkXFJk3SvPB1VOsk20fM8d4c8WGy3xGzPcGtltrZ4VMBqdUpZqDFZX69smyCiRJSEywec+mUVIhfNulM8ViM9q/FTXXN+t8+jbSGysxG9eqtP+D7jZKHXanTWHf5hNCcyB7FHXAJvs+mfOyTD2hWU0XyY59DjAoZW9R0nyVqsZ7aT+id2nOyyFjSM+E5rVeTEiDu0s0m063aqxpDr5HX7+XM7l1QBgq/t6uTTwS3e27QBae1VobxXbHzJQ2hzYepydsqsw/JV5/excfWsgbReqe/rYGOKLHLelSZNaWt3nOB9oN/v81tvghLezZszjVHbtzUvJHkiZdjyk16EGgEDTNJKjXnGLema40JjYXmYMtN9AVVH77qv5UXfReqeHC0/8Wt+o7SZbfraXlMQQUcd5/X4VX65LZNcg4OqZvv6syHFS64r2jDwJDfK3JT1gEl4hLEN786XsBSugy7x5brjwxG8F106ydn9XEqwKLmnZ6ez00zpWEUtSoEsDXCq/gk/664jL69RjJA3EDuKsRn6KclPqd8TdSDwZxHeB4I+Z9xFaKltbn6dF/W25X3NwKf8oddrTHD/5uZDe+mUC+oMiQo81goZhjjynahAA4BwE1ECi5+2ryN+RfmOx0ct5ANzSvjxqWWTb9NVVdqvRPdmc5bX/tQHHSLTZEGwju9cW99tOPqGmxqtoPdu/XyGKckSR7ftP1VHtt1ceTc5HTtjvqMXGmx3Wj1ZdnIg+S0TfIaI3iehXhOt7RPQvieiPiOjbRPSL62/qyeBMIX2SjSyX6ZSumYIrBZdQngIuUnNpwTquykMhqYceqhs4vncU5Quy+zDMwCGbR+zYsVc10edzcYlwEQazRlWfiEIAvwbgpwDcBvANIvoaM7+uFfslAK8z8/9ARDcAfIeI/ndmnp9Kq3X41GM7GCUo+5ieEKPsdLYkpbkwUrSdm9cMkhCnAk3JNivytyJGUVgi/KK+NRJsGXeDE3T0TRoTltkc9LxXLNpI/B8G8CYzv1US+asAPm+VYQA7REQAtgE8AOAJAF8zpBcoSDQK2ST9qliG9Pp57Rq3TbMV5vXv4nl+ngXGAcBIzBEQI1D/l/Uq6Z7nhDwnBCGLmX+KQnCSvvK1bzPtkW5eAhdBaraBHf1sJ4k5L7Qh/rMA3tU+3y7P6fhVAJ8G8D6APwbwd5jt3DEAEX2BiF4jotfu3r27YpP1CrX/dUlsZ6oN2TSKCW68xnWmQmXV57XVvXB3bJfjkI6Qay6yqn0U8+K6Xc41JbFOZ1ZiD1vKZ6VtQR8AcqYq6Ef36zcHAL9PMunFTgmeFALnjlVTu53Xd2lDfLc1aYG/CuAPATwD4AcB/CoR7dZuYv4yM7/MzC/fuHFjyaaqShq89bSOziQ4vigCOYnkqF8RXjrvuqaeox8AEDEoQHH41H1XfEHb8gBm056zuBHJ19agaH8X6bFNdpMlBwcX2S/KAHDSdpzH92hD/NsAntc+P4dCsuv4RQC/zQXeBPBdAH9hPU0UUPl0Qpa8fuG0PpD1v9YWjljepNI+p7df15EEEts59yuwXF4hSSIkiWnOEcN3rffIOWo5//Qv6g2jVYNDy01D1rXJ5lmT6KIMPsuiDfG/AeATRPQSEfUA/CyAr1ll3gHwkwBARLcAfArAW+tsqAjPfBOAfy84slT5VhLY8b8Oq7O3z9xbHBTmxVFa8CnkhTW/JXR1X8/dlyQRmMkZs885VeXdhNfKB2xu7Ck10V72W/OAfJ5z5k0lPdDCqs/MKRF9CcDvAAgBfIWZv01EXyyvvwLgHwD4LSL6YxS945eZ+d4pttuEI/KNA62rVtlj9J4nLwlQyP50VKLK7+7RlYLiGVR8UXkU8MIN1+Xmqq67gn/K0yorTxgUz0ulWH0y30sx6DA4dwwYpb++KP1dW2xJz1yRSZtMQIWztvK38txj5lcBvGqde0X7/30A//16m+aB9IYs8ptOO5p3ne5tp0NMgik+XG6TtCuOjTAHBGnrJb0YC+AhiTLSJQGiflqLtFOw03JVt5f3S4MfBfmC/Hm9jHPH3bbk3yBs+mCzeS67SnLXzlsfc4/TTk7NS02Jy4hYf77dCbwjt8Nzrg2M4JsWEjKdRYg8abVrwTwtyhjkrwo5yK+3T5FfnZLeURe0c2bYPOJLMJb1FjoT5ZDX7Rv6Fuel84zYOZubUxsIPBI9HqbVertCXkpRXdob6r46F+Zgx3xdIS0TbugDgIvwdh4/CuvPrJ7nWwVQRXIPkZ3vd33kP63kHq32a1jmOygnMn16dcpq/2YRX1fZq+U7oZw+YZoHgO1GynDPhQ3LOhqNUc5+bZNdeF48LMgYhLlBfpfzjGTga0N+wBwAJHdgKXmnveee5MZLUTn/d2wMUpFA3Eobp7L6cppKQyvS63/tm1syeiktcgVsDvGlratqGz/It4pRbcKcvG7JLuqsfgRiUIN6XD3Ho9JK2W518tsDQVUmYuRSfr6WUYYq9l6a93u33Grq7VJ9taXJzVbhW2dxPsn1M8RmEN+V7IK1/22Uo6tSO53kB4CgNFIJz7F/cH35iqJSOWvaZRaLXHfOiktIW2e1QgP5pYQbCmGUi9tqR2Wbk6mnm1S/wZM5P28t4dcBzztct9V/M8Jy20icBgc4YhTZbO38cTnq20q1eCzptgNhd5paE+3lQeJKvQ7CHEGYI4yzas3eNQAEvlgDRztcpNddckOrjD4dUO2rnlF9KavCCyTR1oG1kH5Zhp3R5gKbIfGB+pAnaY/6+yKUa8vCS0zJdDwRntWa9DpUnfaut46sOJFvE8yAgZbp86WIPOXEY5NexdzbO+4CQL80AOY5iev7xe69BITlmr7DOCfG6p/SfL4tlpWYXtIvS/gA7u/uOn/KGtTmEB9YTd9pcGxZbi88ArWJqiufZRO+KX2WjTDOjPz5YZyV3nWLMsrgZ+fZD+PcO20IiA3yR6E9QBRtVwOAuBKgOrd+K5N7Pd8mv/QqTrGzn5mTTNtUcL7zpzxt2iziA+av12Td12EZ8yhvkPrAcqm5LMSjRDSiKUNZFGcG2dIsQBTWRXwUZoijDJNZPdjG58YbxhmI2DAS2tl3VHSeKysPUAwA6awh6WaAmnZikP8UUnCviibyt1LvXYbbJrW+LfnPYJ6/GXN8H9oajIlrL53yFm+yYf4eDtIqK20QZ9UBoLZsNizz3wF1CTvoJfBh2C9ymrjIrm+rpf+/srFQQ5tYgWr6o+aovvfadrA+Jdj+/SuH+xJW/w4kHHbjBKxLEdg8ia9jHUssehnHW1XqfRDlTueXpuy2ivSj4RxBoEt6bf/6KENiza3jaFHvsD+vJL9NxDjKKo1BrxNwLw8qSR9q7bGfr8N8Zu6e5wOVyg8Ian/NFqN/PpvVgZWkp6Rpovxfn8KpPmLv88AkS337++uNFLCO17O5xG/68h4nHbOc1gOsAaBpPq8bznS/dd3YRgQM+m5pbqv3Evl1DPtzjCf9qqyvzlD//lEmGvQCYjCTcx1fyr67GPwc6tayPfMcyb8WSDkfvMZBtFP5T/E9bCbx1/QunO/UIwqU62zcT5371QN+skvYG07xeDIA4Cb/ze0jAMD72WobXYRB7gzOkeAnfQlbW7JeKqnNSZpsJedE/taQHMhsSe+8F559GVo89xTew2YSvw2alv5bvEul2uqd3bsEFzJGw3p+0SYpvjecGn8Vnt9+hGGY4O2jK8b5j1+9j7ceXPO2PXTZAkq1PssDQ8W3pb5O+iDKkaeBP6iHitRhPHd8zxMYSs8cEtmY5F2FpRUM3YZQ/g41j0/CuUYsbg7xK3W8qVzDdUY9as+a8JHl26+MWzbpldOLqFpbUGp5P6pHywVg5FbDn99+VP3/4vbDivyDsLj/41fv4+54G4fTQu3XVxCkdtiq/FZpLJyWmXnUdZd7bqtIviQw36N9z0Umvy3R2xgAyoFAjAK17ucoBwnekfX70MrCf1JcfKv+Mp5MLd6RK1RXWaVt0ivE/VS0butz9CDIDcOdjVGvXbZxnfQKL24/FMvuDGbYGcyqz77BJwpzRGGOQbx8AuSwl5mSzCCwx7ytn66SiTgeolu2L5KaD6xneinlX1R126/zlHHxJb5hfHOV8VxTRQI4jX3KqUbyyIt6ddVezfN10uuGtiDIkVtzaUV6e57dLyW4kvoS6RVe2n4AAPhgUstjit3hDON5LN5nLx22QRTlSDUJlTvDcRcndAMnNa3dn/N6/lLwz3DqkXRKE2h67U0GQODU3tHFJz7QvPbieYHOZBzqVnvepg0iEumr+zySdaStyY/iOSJLCziYDfyN0nCrd1D9f3e+AwB4eniADya7uDE6wt3xNgAgyQPRllDzyNPaEocZ4jCrpgs2bPLLaOlI0UaKqVfaYknrXKB/VVfotL4dWtBA/nMc/DaD+ICb/C5jqe6j7ki35XSfbdGXozCveeaNGpxwFD62+6DymhunpkfeXjzB42TYWMfTw2JAuDveRrKEpT4WvAO3+nMcC56BQEH+eRnLrzSdej5CblbjffB1/gYfi7Whbf/S+oa9ExIFdfEvJTOp1e/7/qdkF9kc4gOy44SE2s4zhdceCVtnOTsUu73elBTVPfN8pE/zAFGQ42p/7Cyj8HA+wpVeUe5afAxADqhRsEmvS31b2g/jehtV3crYN09W6xK1zm1/brNu7cNJl7VciTGawFhYwlSXiXOPncL9hdRAUQtiOgcL/2YR34bUcVySfckfqhbV1jBPLnztzTKjeGHMc5F+FM1rUr/WFiugRl8BiHWvO4fkt6cayyCZRyCCM2knuK4BkJQsRLxX3SCck7CK9PdOEZesjwAERfShuGdCS3DoiGB04RSk/uYQv7bcYv11eOqR7jZZ3eP/0XxJKySEQY5ZGkFtFxiFuUH6p0aH3vsV+QdhIZF1qa8wznu1Jb/3xnvGZzUIxL28piVIKv4yUNoNhVxL2iGmI2+S0K6foI0W0Eb6L+OP27aslItAvRf9Wouswhxq2ZyIgYD896w5rPDJWc5bItzVB530kn977bGCNB1ESbEfXXm0ei4x5nlUHR9NC8v9UdbHUSYb3wDgyrCuScyzEGkeINU0gDDIxcMeENp4HHrTgevS3vW7reOn0oOB9D5y0kQWbaePCozFdEDfZ7EhuIsCLuxQZA0YZ4SLT3wbuuRW76xhmU6/l6J8sUmldgCypFfkT7MAzDAOaUea7f7M+DyIUjyYjWrlIsqrA0Al7XW8fvg0JtliiW4nWnj22dJebYY5t9x50zzAQJjb61iG/PkyKipQJ+FpWrHXJRXtdfUldjGqVRXliyNY7Irk3LClDSPXYOjcPOK3bPEyCS8AX3oqv842S6LqsDHQvPR08keWj2cA9s7zJ1lsDAAAKome5gF2+sWAkDnmgUczt8agEIcZ0ixAHGSIgww7lvtwnpFB+krq5wTkVEwFclo439ixrkoynjbWORX2+CUZxSKufTcxlFkYQOyVAW9bgLWtbmzOHB+QfwRGsQuOcpqo5lzq4gKLLLOWMcpl4POQ3nbQAYCDyQDb/ZlBeOf9wg+YM1XLfDb6QYqcCUke4v1x3YFHISwHlUxzYAgpx2QeY+hZeQgEVu6OpsAImMxjTKaLgafK9iMMNEsbriScdH173evjwtehnAA7gWoJX/i2ZKjjkK0pEiy3Xax9OXNzJL7vewdc9/LM3MO1sVGFg2jxwE1eifQK94+3ME0jTNMIAdg49uIJ9uIJtiJzOnCULiSyPSD0gxT9wGzLM6MDNCGkvDoUJpZnn942G1LufQXOqZXvvlnhcsUvBHwSv8Gqr6v1tfMuG4mkFUlOTWvAZkj8BtK3u5lLLaD81GDkTqYRol5W++F8pO9bPvDjNMYoKqTsi6W7rcJOPMVhUvfguzUoSK3W8HXD3qOkbitoA53E0yTGVm8mBuNs9eYYJ+aUQw0Ww0FiSH0AQJQDgmdfTeqv0md9Lqu6RHRJ9yap76rfpVXq5xXps6BxSzTOqd1Ox23e0RpDdC8+8VuSXt8rT5rfB1ZSDXu3Gtu1cuFvXTQgDOv+9zps0ivo5Nfx8eFdYAjkHOCtyXXjmiI9AGyHC+1AJ/4zowNR5ddJrtbemRfnB3GCjANDE2jcMEOD2pGnQkl+21i1FpUfqJMcnv9rPh3a//q+fXb9dn1GHdZycIOkt21F0tZnLLlB29/B5dSzJvJffOJLP7oHMumb17B1v2q7jrDy1GORJC7SK7y4/aAgugMfH96rkV/Cc4OHuD29Un3+np17eG+8X7at8XbRum9/H3158njWM8J1h4MEhzbxgYL8ruQga1m6W6KcK5X1aU01sgCkciyewPpfQRrETsFdudUcn4g+S0TfIaI3iehXHGV+nIj+kIi+TUT/br3NxGL+41nCo5Br66c66V3eZxTkoCBHEBWHbzspIq6uX90ee0n/Yzfewo/deAs/sfcneLF3Dy/27jnLfnxYXNOlvY29aIIb/SPc6B9hGJohvre2/E5C9nfPOKiRfpa55YD+vWt1u/bNq1WiHacBXzN8z5Su6Vuh6W22Q2tbkF0fEJS0b1x1UhrGeaXeIqIQwK8B+CkAtwF8g4i+xsyva2X2Afw6gM8y8ztEdPNUWgs0qloVAkY0MCWc7pBDDW6sVO50Ewvaws5wVklG5eOucDzr4cduvOWs98XePbw9L6R7boUOfnx4D4/TIkBnO5riKF3YAPaiCQDgZnzg1A5ubR3io+Migk8nei8qluiAwt1XqfnLqPgGdAc2pbZKbqVqz7xl5uC+ub0Lrj0UVxHAOpHt6WDIgLU3AmlzfN98Xt7tWCur5McZ+ey3UfV/GMCbzPwWABDRVwF8HsDrWpmfA/DbzPwOADDznXU3FIB/RBcGhDwNjDlXEOZV1hxmauWM0isluk0SO3WVws+89IcAgIRl1fdHB+/jRwfvAwC+MX3KuPad2dPY1px0tqMpwiV7ryJ/z5GIMxR3upCl/VZ/joPJAGGQL4J3fM3RyG9umum4r0kKt/nqbaS8qmdN3p3LbowCFH0RwIXJQtRGR3sWwLva59vlOR2fBHCFiP4tEf0BEf2CVBERfYGIXiOi1+7edc95zZtgukI2Qe1HJ0hqe38413bUClujhWFNqblts9fElCGmRRt0wit8ZvCh8flT/Q9q9QSUV4fCdjSvjk/t1sfYW1uHxqCkpD1gDmC+Kc3BZICDMvmnIv18FhlGUtFIFbB7F51l0XAPxwwOF4e3HpdKvqrSs4ThMm/IaWC8xzNaYG8j8V0LHHY9fwnATwIYAviPRPR1Zn7DuIn5ywC+DAAvv/xy+yHTYeCrEhmGpqVTJ72S+jbpq7Jlh7Glv6Ti+4gCAN+7V5/D/8TO6/iLvSNMW3zbd5NrGAVzjPO6F9/j1L2Ud2t4iLvT7erzw1kxXVDkTzkwvAUV+YkYh/NiubCnue0qwrsQxDlyx7y+dWReW0iSP5A93li5X88D8bo9IFWrQBCe0QKcURFcY59D+41M5METp67ytyH+bQDPa5+fA/C+UOYeMx8DOCai3wXwAwDewDqhVlak0VsjpbRRZNPvGoSMOMqqHWxtjHqJGOGWlNZsifQ/tv2d6v9B6UQwtdYNPzP4sKby6+T3ER4AHjUk7VAOQTb5bdVe+fhPHem7jDpt0kvqvYLtgQY0S1n7+rIrOr28dOByo5aZqcXUQlw5sbZmq077ArxKdV8kfVVGVeRv06poo1h8A8AniOglIuoB+FkAX7PK/J8A/hsiiohoBOBHAPzJeptawLvfnWdK4OsGYZBXOfMCoX5fko04zAzS/+lxYdfUSa9jh6La8VeG9/BuYqbLHgVz3JnvVtF5UoSeRPqDed85j0/L3q6T3vbak7byCsMc06Me8iQwSV/66TuhX6u5oArlW1r8RWluZ8PxqP7OdGzS8yNeHC743oO6ph08D2XSs3Cc0rbZjRKfmVMi+hKA3wEQAvgKM3+biL5YXn+Fmf+EiP4NgG+iGKN+k5m/tfbWAhBTJFnvMJ3GNYs+YDlfeaR7EDCiMGuVSkvKfusiPQAkyBELo9MPDt7B/zf5mHHuhf4DvDO7Wn0+yvp4lAyR5oV0Vok6FQ7m5uAgxQMczvuVWq+iC+303rp9wLsBx0mdSWyPuDNAUw5GAEWb4nz575YGhU8DsJQBj3LI77Ja0gOWThrSgFYOPMz8KoBXrXOvWJ//EYB/tJZWSdB/sCoSRx4JgyirWfSN62XHDgJ2kl9PWa06f5KFlbq/FbtTZb+XXMGzcX1AAIBrwRYO8olx7o/msrr+59NrRtCOkvBRkCHNw5q6npQDgm7MkzDPQmNO78J0HiOZlvn2pGQbgJv8PmkPmBb3ZfpyP6sXd2iBQY+RT7R3NCzv1dpW80HQPTzL5Tsv7P6TFuHbLo2DMzLemUM5c2NN0v/ie+6t6pQBzZoaKaLX37JN/ir3XBqhV0bZKQmo5shS7jqFW4PCkea9ZOFh915yBT+/s/DV3w2GFfl10v/Q8M8rqf/nU/9OOTbuTxe2ADUAKEghw2pOP5kt5vTDfoLD44Hzverk1wlTdV4VJdkm3559/QTWbGmNvLrWa/DXUESfBybpFdqQX0CVZlylfZPMUr6mud7ZpfHVtyQCZVTM8QxvKBJz4iv4Nrkorhf32s44OnTDmEqeqUMRXsJ+OMa/GheW8s/0C3vAgEKnpNehQnV9Rrw7021sxXMclwE2EtEB1FJ36aQ3PpevMhHccylk8HFUHxsUOTLTYYdy8ttlFo2rk7+MutSfXdv1qBzUbfLXJK6R6UYYnGzS620JuVopyqeR6Q8gDAp6jsJKS2qj+vvW+C9lWK7LGKQXEZaRlAuuhIC4dkjQs+ToUEkwPpzs4N5sy7j2R+MXq//3w3pqrI+yEG+nwJxD4wAKqd9G2kelOn9HW8YD3KRfJ1rNk4EqpxzlVN9Oq6oM3t9X/FmWIAG5ltWqpeC8OBw7IwdxZvqE2ANKG0cefY6un3a9R7vopd80U0l/hwRRkp9zQhDn3u2fJaR5UO1tl+WBSHgbSvLfny9U7ZhykfBAQXoX5hziPxx/shZ+q8frB2Bsx4V3XxvSu+b98yRCGLCRtSe1NuNgYUmK58U5Y7NISeoJRDek/xJTVZcZwY5rr3YyHha/Yao28HT0AXtQoDivpjCuwK4gzpALuwi3grjtTvl3/YZ7LzaD+Pa7SgoHDSmhge2Dz0xIs7C2D70NtZw3SyP0o1R0x7Whq/sfTnbx/OhR9fmbkxeMst8/fKd2/83wCHeygrz/4fiTjc9TOCrj+O202z4keVh9pza581Nh11tF+uqztFMsiumYy+1LheqKnnYeQ5/OXbJIaYdcK6idkKIow+RoMXhG/RSZY1ffeCsBESNdhtxhbhkzaaHu21/TOYqVf50xD1ir5N8M4reEPhDkSeDsEE1Q5PfBnuM3EVANBOO+ueT2nenTtbJxkFWSWpf2QDFP19fe1RTl0XRo5PWXEoGq1QnV1oC4kvpK2kuE9yIjOe6+yRlG2WoU9E0rHFqdTfjISn1lr9BEWrzCcHtmbP0d9rIa+ast08jcDj1NwupzmoSa1PfPT4j8rr1OZdQ1119jIo7NmePrH8N2I2CeBJWrbpqFSD1qtg17jbxWt5Y+WxHp4XxUHW9P6vP0jAOD6Pr/tiONb0kuB+FIS845Seredmpn3DbfgQjIsgCZz9usDEypZbS2jarLqPDOTLNWJU0ONOo27b7IEaSkI9T2Rqztk6g1TR8EojgrP7f7oqLRmQrjtDe+wBlbsB6Hns2R+C41sGlDTQtK7c9y8m4pbVvAATOyzZcYEygGgReH9wGYiS8VJEmvHGnePlosBeZM2O9Pa2UBk/D9MBUj7FTiz2kq/9R2nL64Xu9bzmrwzVe++66dZyrV3/VbeAifzsKa1A8CdpI+jjND6q+K+bjZrRlAtYzn2z/Pua2Wwim57G4O8QGAUexbBkBt3aRvj8VZYBhsgoiRJSHCcsRWLz/lsFy/X1SdpKGRUtrOSiu5wTaR/1uHzwAAtiJzmfCd6VW8MHgg3SImvnwk7K7r9agrEYcZ5iXh7cy/SRq6E5MoSTQTOmPEC7Krv45BWQ/YaRoA9EzJxv8eqISoesShHlE5m9e7t03+wda80HwsbYeo6FcSYdVKkTfqzp6me3wNAABROeAq37QcpxrCu1HEZ33OTmXyTCGLqY0sCatoKWXlz3Oq+eVPk1hMT+XyfQfc5I80Vf047Rnkfzwf4o/nz2IU1f0GfKmzjTaVNgbXAKA8DHtRWpG/p5F/EKU4mDTn23fClvQMBAkhjxnhlJALPYtygOZKwguDgE14xyBhZ0COy1DpnmUD6PeK8/YAoMivR2Cq/pFngZmF2SJslgWVViQOAB6eNpHfcJAKSkXWzhu4JmwO8WsEF8qUOfGkJbw8C2qhki7yA4v8dMMoqe1Oo2O3N0NUxsuneWgQXoci/+MGp50mLUJdU3YFNQBcHxTLh/PSKKh78gEm6au2Dwvp6BoApL3daV5/8UGi+flPS8Nh+bg8qlv+q3DYVNh8UncEslR2TkL0dszU5Ao26XWMBoJj1gAYT+vhz23DaY171ADQYhoR9DLkmlFR2lob0GavuufgGiP2NoP4PtILYZE+8ocNRh/l5TfPQvS5lCTWUuAoSmrkzDlwkl7h2NotZ5z2DKl/+3ivrEsmv2tAUKTXcU07Zw8CNq5sTfDhrD5vpYgNxyiJ9E0Iy+lCHjt8L3xVtkxe6SO9muLoOwkPohTTNMJoMBfJb0NJap/xk0JGUCZeUcSWnIc4J3MZWkjH5U3y4VnyXAabQfwStSUg42JxXs37mQmc1bPspEloTAeyrFD9sixAX4jGO056RkCOlCpbwrTc8spe5utZg4MivyK9gn3f9cEYD2YLbSEgxtW+Geyj6ldSX+3Hd6Us93Bmahu6r8L21hRHx3Vbgk1+G+FYVnOVFqCIrWsFWX81q3Q8SowkIgrjSQ+jYV2qSzsa2edGgzmOGqY8ytU+1IlcCm1pMAjK/RikZCU1dT8oF/v1fJCu17NGJ5/NIL56UQHLKr4Hzn3dNSRJYeyblUkopAGg1qRSKkeGpT9AQHlFegnvH+8aS3XTtCgrORgNowTXB8dV3TbRnxk+WtQ72a/+lzbgBBYDwL18S7y+DPkpJQRzORhHJ7mEcErIe1yVFbWBXl4EzjiCbJgJg35Sm78PohTDeI5JUpfkrq3Ntsspjz0ABFo6dSmS06UBKMHiy1RkVtRQJuTG1ZNlsRnEV8gJnAuOHQGcobpS5lP7nD3PVwOACs1VUt+W9jnXPdR8pL83MQmnSG9jKGgVAeVGVt6nBo+N6/og8OFUNhD6bAcKNvl1tbMprVYorQLY0Jqg8pE6yQ+I5Nd38725vwiO0v00hqWWNkl6eHq72J3ozth0cQaAXW2pNCTGoUf62+Tv9VLMrYFnpdz6VtCPmGiEsVKUoAubRfwW37tt9lNFfinjDgBc3z0yPh8nPWx7YvCboJM+yUOnNV4ivYIivyJ9kyHQvHdR7uboCNMswkG5TGiv8W9vTfHozo5ZAblddFsRHjBIH8wJ2VAbfBvIH+/PvNt3A4XWZDtp6fsN3BwtflNpiRQoUqcDcA4AivxK7e+VKwfzeSSSPojrG2hSqE0D1DyiuigY+9T9+nLnCbE5xFcCPQ3AunHE1jbTYh2kzQDAgmFQ+ew/Hhfz4d1ybf+qYEADgGkWGaQal1JcaQcBsShpdBxOi052bfsYsyzyeg3akl5HWLbjeq+YHtybb2ErmmMrmuP+rK7e7/anSPMA0zSqefnpASs6Kot87iZ8UPIzj1FtbRZaxvj5HtfmrGrqkO8Lc3OB9Hujuo1DnzLt9oqHSgPkfn/qJD9QDAA+6W+j10txoxQW793bx0DLADUZC6sH5YBgDAptrfWXwrjXRnvU85xXRj61Y0lelTGqLctn5RpsGC3y7rngylwjdSw1ALz/eK8aPHQosivo9+seePZ8/SAdYDeq1xdaz5/kceU7cJz2qnr0qUjq0DqUpVuFqtoDQEV4KYZeQ+AT0L5xeRwCI/M9T6YxhlY6tYhyw/D6eF4QWRHeh+1o5iU+UJD/YFwvo7TEvVHxO+yU0wX1uw2sdg5Hc5H8gDL2wZDyhgEw19b1c9RWAVbFxSe+PgpW67+l1BffgXmSM21TPN9jckKShiL5XdJeJ07OhKnlMnv3uJD0erpq3RglpfDWcUWw2gPAUBsMVO7+iZaSe5K77QyDMME0ixtJr0NJf1HCq1crVBdOURisBdeF3gFhvissuZZGP5380VZ9inVtq77V2F5viv3ynR0lsrTe1rYof2brAO8f+x2m9rfM32CoDTSSAREoNLf7RwsNK80CxP1UTGwClESXDKI2ydfoxbcZQToKOVol5bDh9Zay9jBP0hCJFpdu55j3OfPoUKRvAyXtHx7X19sfWCr6bjTFRDAgDoOiQzal25aMj1eHYzGaT4crUUWFvDj694H4oCR9iVAYv/IQiI5X68gS6ff7k4r0vSDF1f4xduNpdQAm6RU+vf8RPr3/EV7YfgQA2BtMjKMJO444imvb9TbG/cUUhlWmIld/tu0Cl9aqr39vKZkDC+W0C7olX6n5PgusIn8cZbUQ3F7plpZaG19IZI+jrKrLXnqy3UaBgvxXtsaVtB+nMcbpPp7beoRpFmOaxbjZLwxW+k49QEH+cRBjLvjL2oSXJH6aBejFqRyvP2ke8PoPF/+X4xCU8hFOFpI/d1RVSXuFcYjoxoJ88zTE03sHRpHjpIdnt912D4XtaIadeIrDRNuPMDaJ/emrH9Vcpnf6UxyWU4KhZdwdxvNquVVB1zT2rk7xZ/evGf4SYelvoOcHAIqQck4D2YvvFAJ1NoP4UgCIy5nHmndK5HYRXiLi914zN8pwrZEfzPvYG0zweCpLXClgBKgvsalBYmwt9d0+3q862Z1ZYXH/3q3FNmT3k4VmoKfLlox6Nh6MLffeOMUEWscsSa9b9cPx4keJ5JkQgnThugsAaiPg8a3FueiYkG65B2B7ByQ7H8H37t0TNSBgkY0YAHZKqa/+6jjWch48MzpojJe4OVyQPaK82q/Axt3xNnaHs+ViIlQsf2jZpVSfvTRzfI+Lokh+ZQdQVuoTeDvFUYZxGlcWepv0o2iOcdqr5bOX6nERP0nr4u/O4x08f/WhULrAuHT9feNowaBr/cKirDq6iwwKvTDzTlt4HoihohwA8SFVgslFehu61NIt/ByV8/19RmCps7kjSEcNanu9gsTK5iF9Z538ADDSfsOx4x3Z5FdSXye88QwP+YF6TMRsFjc795Q7Ddek/Zoi9i4+8QHTgNS0k44Ne50URYKOqJ9Vbp+2B5Zu4Hs8HQKldmi72waUN5Lehu395fIjUNOGG1umP8E4bfYtP84WZeygnqKOosNL22RPNOMeh2yQPz7UnFc0jdueWcTlYKCa4VJVuan37SSYHfXQ364b9xTpdVwtH2zHTGyHM8P4qTAKEyf5bS/Jp4bF9EqaRgF18s+F+Ywt/Vt79p0CNoP4Cmo5Q5MElBEwaJgEWeS3kzeEYV4lb/Al56g1h+tptq+NjnF/rDnrKInuyrPQ4IRz93gbN7aOME56GHuce+7PtnGtf2SQXp/X6/n5ZprDzv3jLac2AiwSRQzed2sQQQoMHjAskwPCKZBs1794/xEw23dWZzwXQI38EumH3rXDhfHTHgBGYSJuRR4HWWUP6AX+bEwKyn1bIr1CkkTGYB+UfTE5Ft5vwMX0KqW1x+VvllXfAV8QSQWhiCTx9Kyzep79faGz5Uw1N97H0yGiIK+tDvjgy9cXEFcDyb3JVs3tF1AGwBjvHl/BOO1Vh1RXzlTLJ6ji1r3wjIeDB/LF/mP3gOyT9pLtZnZUfB+b9MMg8ZJ+FFgGOetzRBm2LGu/Wg6V7AGjaIZRNMOV3jGu9EzVXzle9YIMvSDDjZGprU2SCKOh+awsCZBMouZIxFVcgX3VrbW2s4A173NldJFSQ4Pq0h5ozs9mG+z0PHuAHLGnZ+/R017r0YJ25KDCoJ94c/2rAUARXoJ+Xu0BoFvylyF/eFwOYMLSUxPpt9/L0H/E6D/yd9zoiJz73Ksl1/k4rr6z77srKNLbA+u1+BDX4kNEtooiwJb2qSXN+2FaHRJujI4wSSJMhJWSzFbz9U1iJAm/RvJvjqqv3kPLtExAQSyVbluFbSo163hczLVcUn/XGpmbsuiOogQfHJnW4GEvwUTYdjoIGbllOFMRgaNeYqjiPihPNVsKKkIcJT0MPO6/dl7Bfi/FbB5h9ri/ILuE8vXbpOcQoKwu6ePjHMlWYJA/mgCTa4vn57tA/Ljco3BvsTuOjpGQhOP1x7fw3NbjmjejLelHls9wYs3Vt6IZvnt8vVY/UCd/X/scUYaU3QZVOxRaoUZ6hYAbts9eD/k3g/gOzunS3t5GS8+gKsVqq9xsagBowqOSZLrKP89DYz4XBrkYfGPn6O/FMhnb7M4r1fl4PqjIb0vBaRYZ5Lc7Ys6EmeVR5iV9ib23ZGnZtOYczkqDap8wvM8G+RXigwDz/Qw4iEG75jv58GAXT+0uLIvPbRVr+Adp6a4bTXE9NlXswNGoUTCvth9/e+zevWieRxX57WmBjd3eFAdz0+lLRXkqN20n6QEAZT/OaLENV8QLB541zfU3T9VvgdBOlazBtqo/d+0hbu0WFtut/rw6bCgnjkfzAR7NBzhKe6IRJ2Myjrbz/GVILw0uj+cDfHS8I5QuyH9/OsKd8TYSYQmv3zcHIo9tqoI0Cxl9MMXg/hyD++b76z+YY/vdaUV6HcP72jlCPejqoC5FPzyQ19n7QYqZJclHwRwDamec82GeRx6LftZK2us7MMuwvPVsjq/RwLcZEt8IW0QxGvbzegQeA6E1h2eHT/qVobkAfWv3EEezhfRXqZpix4460saZzuZbqcCiMDdcZG3S96O0UvdVos95GiLjoHF764eTwhlHfb9jwZ9c3+7bRnq36Kx5CLgyiY0+Kv7qEaSjD0xVe3B/jvCo6OjpTvFeh3eKz5ObppY1vM9I9tydmgXJr9C31PAbPXnz0gGlmHKEo1IzUPdthzN8u8yGvCq2o7kx4FwfFBqHpOYfHQ8c6bXk71+R/zys+kT0WSL6DhG9SUS/4in3GSLKiOhn1tdE+yHF4Qy7NXZ7AG5ePahJcpv0Ctv9+oictEhjDRQ72QBmfr5jj8dW02YXgJzdd56FNcebyTyu2RIeTkYi6RWSLMT79/cqNd+W+gAQzAiDDwPEj4p3MPpoQXodO28uVG9KizYr0ksY3pkhmrB5HAOR7B8DoC7574+3nIPguFzStOf5uuSf5VF1tEFQJlR97IiFuNYz+9Q0i/F86f8vgUKujrXE2S6Jxm9NRCGAXwPwUwBuA/gGEX2NmV8Xyv1DAL9zGg21h6imrbGleT0AfM/V+5W1/IHgXqvIryQnsFCtdyweK6n/0VGhYvu23VJSX8/4enP7CE8ND2vbZAXEeOfxvrMuHZLxcFko8vO7I+gesr2HRYcclF7LthYwuLsgd3hsvm8O3Z35+FnPgFiSf35Vvq6cq57eKQabd8f7xp6FQPH+pnmMaR7jallhwu75y9QIg7as+MopR/DM6wepoRFe642NzVNtHE77YnqzsJ8im4WQBgA1z6c159hvM9z9MIA3mfktACCirwL4PIDXrXJ/G8A/B/CZtbUOEOd9Cmr93h4AdNIfTfvYLudWT22ZauBVLfrKNopJLq13j7dredv0Mvqeey5pf1XTNpQ32HY0q8ivBqWtXhlLP5el9v2jLYx67oxAR9MegB6IgCtbpjRSEYf9floz7Cl/fEV6G2r+r5NeQjDW2qaNmBzWCTTfNc/N9+Q6XfkSdPLXYh88hP9oXrcV6IOANJXbjSe16YUOm/zPbz/Cu0f7zvKAnpefjfx7tXyRayR/Gz32WQDvap9vl+cqENGzAP46gFd8FRHRF4joNSJ67e7du76iJhpWMHQHHpekXwW2KpnnhLEmYSVf91kaOUmvjIgStqOZuG6/1ZtXg0D1jHJN+GgqP+fIiql/eDzCw+MRPrq/WwszBoD0wQD87qKz+lJeDx8wrn9rCo4CcORINjmZA1puxN5HxfdWpB99uJivu0gfHVgbWl6d4uhY/r62xFfYCadIOBTJr5P+xuCodl1lJpKgVhBcsCMhnx4dOlOt+XbjOU3jXhviS0+ze+g/BvDLzOy1PDHzl5n5ZWZ++caNGy2bqD3RMwBwSs68xEfTPrZ7M9FvXYLucivNI8cO9XoyizGZxciywAi11pueZCGSLKyk/TBMqqMJsyQynIHUd3twsFX+36uRvmpbmQFmMu4Z2WDSB/VO3L9vvp+kjDYeluv2x89oEtwiP03kgTdIGeGs/i7VtHvwwC3pffCRXkfCIQ6z4rsepEOnt980iwypvwz5A+QIkGMnmlWHESDUSzDqJbh5pfjta6QXnM6ITsdlt42qfxvA89rn5wC8b5V5GcBXqRiirgP4aSJKmflfrKORBixbiJ0gQvmdD625NFCEzrZJy2RDUvvH87gy0E20zShSx/Jd31q7b0N083kyoVWW1wcHWwAT+sKuMVLap8m4h/xx/XzvEVXZb3XHtqHDQw8AkqtD9N5rjokHgHCWYXqjj3hc1JeMCvJPhTl9dBAi3c0QXV2Q+Oi4j+2t4jeUSP/RbBffO7pTfY4pqwgPmFMBm/x2BqXqfBphEKXYteL3D9IBbvRMbSHXZOlDz3x/MEgwncYL8jsSoXBGoFPYP68N8b8B4BNE9BKA9wD8LICfMxrH/JL6n4h+C8D/dSqkVwjKqLGBW8GYTHsG+SXoEVWjKKnm+VKSCrX+Py2lve55J6XQkrbsUtjpzfDBZBdPDxfW8EkWYxTOK4u0glotUMt/+mqAndoZAGalxK/2fi+LB4Krso5sK8fwvcWgpTjh3emmRHiUINsbIXxsrZZMZ8DAVM+nN9yGvX6xuTBm5SCQ7pV70j8YIL5WkJ+5CKCRSL8VzXCzd4BMkwxSAI6E/V5B6g8ndV+IfStST83xY8oxyXoYhs3Ty+e2H+P2UaHS6NMttf1WlhNgaQDVkp/6Omt012/8WZk5BfAlFNb6PwHwfzDzt4noi0T0xfU1xYOQzUMh96tAE0HtVWG0kbVUdn8ywiSJxb3mVV1Th4rftPWykvZ2wosPJu6ED4+mw4r0OtQAMD7sO0M6s3m9PfksRD5bnJekPYcF4XVBSDmw864nUcaRoLkcHQPzkgzTGTCdIXrvvrMOSdor0tu4slMMLv/57rPYimbGISFbcqnsqeFhNQ3T03kB7pDciTVY7wrBPQqSjUX6vfR1/ip+wWPoXhatFjGZ+VUAr1rnREMeM/9PJ2+WBpdFP+TFJW2nHRtKzXchohwfWemvd/vTKuc8AByWyy8qfZftZw8A81lc35TTI/UVPpjs1lTMBxO3igjAcP5R5A/KKY/UiSpMQuSTEMH+HOGxlYOgjLXPBma+PAXdLLL1/hxcWp7SnR6iQ03iHan1uAToLQZK3hmJ0l4iff8BkFrz/f3tuu/FO+OreG740BvWPArmGGfL5Uy40T/C1Z6Z7QhY+OzP8siw7E/yuIroU2UkNV/Ka+j9vU4Rm+Gya0l1NQLa3k+9kX/erAxrtrSXoHZYORS2lKqalQXVYX/WzyujnG6YezwfVIcE27+/+g6zSMzWmieBuxNNQjNn3od+q7Q2JUb/sRxTryPdKSXekeyBwzsFCbb+7BGCWY5gtvhuV9+oS/ajT6TgiI3jeGZK1Rf2HgEAbk+uiMbaUTCvHHh2wml1PErcg+qN/hFu9I+QafWp/IY2lOOPlNH47nQxWKQcVAcAYwm29ntFcip4A2tS9zfDZRfwSnV1bj6O0dsyyX/naBtXhpZRRkuIuBu7jX0S6TmnitCr4IPHu3gQp9gXNoMATGk/m0dVdiAK2Jme2QlPsI1rVxxg4UTTE5Slrffl+Wy600MgEf9F2R1Wkf/46ai1+Dme9cQ4CqBYqXmqXxgYbY89oznDexhnfdydm3P5G323Znizf1hJ/lh7aW28/qS0XIr8x7D6F58N6YFNIn4gx2rbA0EyjRAPTAv6w8mwRn6F97S86mqtVc3DlEHPlR7rJHg0HtbI/2Aycm7b3ER6Lu0M+hYCeg47KW+BIr9S82sus9oKSrJNiI/c76H32hve9kk4frpcgbnDmNzUknc+Do05/tb1RcMk8ivCN0GP0lM+/XfnO3hxeA8fzvfFe57uPcZhNsD2aIaPZqZN5jDtV5uWKOjSXsftgz1MLfsREde31+plYElzW3MX3Azie4jHaWDuN65hqHUQm/zO9VnN+FItudhTDX3nngakJSEjwfL/qNyma380wd2Dempu337sRnsajIvAIi97OAlqKbJcPvLzHUCPeRk8ypCNimf175mGAM4yUFhvB31Q+Pvy00Ws+/D2ASbP7VaktzG9Lr9XtWvNx/cKI+FOGX///KCelHRA5a5BHCNEjlEww5TrKvlOOMXOcIoMhKd6jwCgGgCe7rUbTNRc/kpvjO8eOvyMgRrpgUXOvRr5o3wh+ZmKFaw1W/Yv/hxfI30t66t6N9YySDJddKo0C6vj7tE27h5t445jswvXvFqHmpdxRtWh4NuSO51HSOeRsZ4/ncWYzmK8f3+vkeTSwMFJ6CQ9zYKahhROimfojmyRlu9+GcyuC9OgLANngjU+DEB3HlQfJdIP77BB+ujxopGK9MtCDQDqf/VZzfclPNV7hGf79cEEAG71F8uvh1Z8xcP5CPuOzTWk1ZnxZKHZGQlH7EGAYQq+NVn2Lz7xXbDjti3yz8dxbedUABiUIbAH1r5pkyQWd7LR4TKeca4lTWBhANCMRcfHg4rwgJkfQCe/y8VTIZ+FbXYGq6BIXzUpLEhfPW8Ib4fqCWq+In/4wQPjfHZ4iHw8Rj4eA5pvPt15gAc/dBVBhuqQwGGZzSfO8cwLi2VAJe0B4FDznMvbOBuUCMGNS3y6jUB3/gEK8uuk17Pz3BXyIVbltAFfJ72CnW3IqVE2eLC2xcVX9e11evVCpC9PZq6y6aSHgcd33ya/Dzrpjfh6pXVoO/UUn8vLgbuT2VMIoCS/Oq1f1oWCJuV18lPZ/4NxqZWU78omvRfWM4PE38sU6YPrV5G+c7t2Pbv/AOE1two8flpYrwaQXykGaGXN3+rPcWdSaGo3h4UhTlfzcw6cmXYAILGyeyry2w4+PsMgALx57Hc1T/IQsTai6dLeuUsRrOma2lTD7h9n6cBz0UC+fGSCRXQqjK4KbbLjZLMQ2SyCtLWRM4qh4QfK0lAkvbPO8nyehBXpSdhLno7CivTAYmqkdrSxjyboOfVSYSl89F/uAnFcHDkjeu5ZRM89WyuX3S8Gh/Tjz2DvDfcOHIN7nlBebTnvzmS7GgR0zPIYM2F5zSa90bbyhevLfzp2wilmeYT7yRbuJ1vODMaStG9S8RU4C0x1nqne3y6lcU+hzdxGSMYpSX5FetcOuYDgEqv/GILBUWkberYdSSLnGQFZiLBht1wdmW/AU8/Sc7OpZy6x2aJy3pF2x9F9ZNJBgOFHUyRP7yH+oG4Ei557FtkHHwIAwpuFhExfvFUrd/ScI/T3ymJufvBgC7tXj3Frp76eLkXdbYdTHOaFJrcTTL2k1+t5XBo69sIJjkr1/l4i24IA4I2Hxff65JW7+LOH9Xx910bHcpQmm8FkYjZolMt6etGQRcG2KjZH4i/znTNCaAXvTCc9pFmAo0kfR9YPkqRhtWafJCHm88givUDyJCgOaZspR/Qf52aixSwJqyOZRc7vWJFeyyxkg+ZCZFdJ+jys59Fz2LYQzlHlvE+HVD0WAMY3QoxvlBllbw0wueWeKtH3f6oiPQBEby9S9+y9MXaS3if1FVSq8D89vuktd5gPanN0G/bg8aeTYoB6Y2wOVNtllt7bh3u4fbhwK/zmR7KfwnsP9+vPUn2q7B/ZJBIluXMtf42xOpsh8W1jfhqAHUt4bTFLFhZ2ZWizfe6TieP1WI+WlvcU+U3p3/zLJdOo0gSMfIE64cuOQQGAqSMmXpD0ebjww8965h52gLl0p8jv8YYFABz84C3s/v67yG/UY2r5+n5Rx71Hxvnxs0Nc+U6Oh5+qt/34pQSw2v74zjZu7RwawVPPbBUWdkX+T2zdMe7RVX5Fft2SL2kLd8oY/Q/K+OAHiam+/+e75jRGxW74tEYbmTBA6+S3SW8pCLVdoVbFxSe+vm+eBsoJ7Fnfz6YhQk/0HlCQXyIjEQtSu/Rmce0D51jbz9Owyu2v6rQNgQrzcQwKGFkiTAO4Hr3lgku99+0yFSRAOgBUenr19Wd7VN2nQmltSKTXwdf3QfceYfysOecNx0Ai+7sYiHfm+O6da3j+urzMBhQDQM6E66WPveTUowaAh6m5enOzd1CR3gdlmJPSoydpiNkswvbWDIdHxXP0qM3pUfOeh95VGr2fXqrlPP2lqP3tczNWmXMyjjxdHGtvgwOcBcZRtMv/mlV7dWSW9pHNw/pGIvOg+AVbGBQl0rtiV1x5SpJRcWG+YxZ4+Bdl4utJOpT0Vzh4vsxirGkZ86uM8Mgti969dwXv3rtSO2/vbAQAH8728OHMbNdBNsSB4LTQhvS///6Li3ZalvnZLKpSmOlZgpIkRJKE7Ujfdv6eo1U/bMLFl/g6cohDFeX17bLttEV5SpincRXIozzq7CwocZltVlLzRe0gWORJYw/7OA9EvVmvU9ICgPYRXEFS1rWEISjrA/r+E+mgPgXQocjfOywTaWwVn48/VojurT+Xg1oOP7nfqj0AEB5FyLbdyw7v3ruCZ7YOGjMpAYVLblKu829H7qW6e7NtXNf89XeiKQ7TAf7LI9mOoMhv5yy04UutpWCQntgUZoTFb3kpffUVcgAniGRMG9xblU88p0H9RUv9TE+O6Mr8W53StZPmtnoJr88VQwCeWY0RXy80LxsAaiWLsnL+r3Ekj+saw+FzxfMHD80KpQHg+ON7CFJGHslEjQ+B4xctL8MG8r99uJD8rjTWKnJO7ZF3VC7D2dmPHpUps+/NTCu+zwUXaE96OY9+Oyl/GvN7YJOIr7+jlGqbZxpF1QCZhCDPkplrNNYTVhj1puTepFP5FdnkdzXT/tFVqLHdptxangzZDLEtETjn9YS8V6zbS3tEqnN5b0F+oE7+6rxmJKcMmF5RGoBZTg0Auk+NIr9S83VsvU018gPF/N74PoJWpLLYPr/9qIqik3a6VXgw36ri7R858uQ3kR7w51toI+lboY3wWQEXn/iuL5pTzY9fjN6zb0uCKmmFFOBjkF4YYX3kV+C0DKzoOcR6Zptqi3O1qYT6nGkqPFCT7pRR3foLTfVX7Qpl8leP6wGh7sZbTk3DOTDflQcCF8Y3dUci81o2Kuou42IweUqug6YhUBr/bMLfebSDm/vmaPPB8S5ulVtTq73tgWIQSLT1zGkW4/3JvjFV2I6nVWLMJtKPtf0WVYi2PgBQ01JIUapwwsrIrf2tPyi0wsUnvgRFSJ0oZcew5/a61FfGtsxQz7WReR44d+LVA4R08jMX18TxaR4Y5Od5AFKfa44c5RKdNKgoaa6W7hwqn+rHxHXSV2U08lMOrwRRAW1pDCB3awHpEHjwqYI0QVqX/jom2pTZEQkLAMjL95AcFqNP0M9q4dYKTYlVDpMBBkJy04C4Iv9RMqhy4tnp17b7Mxwf9bG1PTNIb7RXGABkSHai4m81ADQR/lIs55WraF4QAyUBpSU+noetw2h1gnPIdZIp77w5gV0SXX+2qq+8j+eLeV+tDKypgm81omxXLWIRQKD75ut2oxZLzdnIv+ynawGpI6pPz3GhDwI66cN5UZfkX5MLg18+CyE1y84+rCS+jmK7cH92JkV6CWpPxeOj4i95tPhqANBCa5mpnOf71X8KCvJzWBr4asFe3tuXwsUnPrD4whIPWqlV7nV2oPDCk6qvSOUwwlTecp5BxTUFKToFybkEJsLSnVh5cxH17iTSS1l4mIA88vvyc+S1JRpQg0DevKKFYEZItzQNyfLDyGehkS34o4e72Nky5/J/+vA6ro3MBAO9IMO96RauD45rm134CK9gb+RRSWaXg53SPi2jHoW5SX7bWcfuCupyyH4hsAI2Zx1fRzUQCIEq9jxZ+yhaVr17lVvPc6HBOitJZaB0tJmGxVRBOwDUMwintFivZ+1+6znSQOMKfwVgpM/WV8fsrFIcLQ71uS2S7VK69xYaAyBLe8rrg5GEyJMu/LEQHAMA96ZbOPJsJKowjBfagWv3HgCtJLDtCEYtNksFYGqueiKONWEzJL4E10vnhZQVNyLQjG6NpNcJK2kd+jTENsDpxeyxyB4IdEOlZODLUf++2mfS7B1qZxoO2Zq2uFV9Vk5AFvIIRaCI4z5pmc/WFiRyZw7uZUPNQJbLP3E+CxGO0ip+4fB4YEh9td/d4+kQe4N6urUHU3MXYbUdubThqU76IMprlnrOCUHEYgKWKnFKnGE+sVJuhXlhW7KnpXamp8Dqw5duOU8yWavzwOKaTgbfKDlbqPa2xKse45LiiuxsfVbQpXCAhdExIVO/Enu1NdDo9Vo/em3wCKSRyWq6z6ovdaoWUka5xNsDgC82Ro+czbXwYnvFJBiHyEfCZiVJiKAkVpYGePR4VBl1r+0v5vhK8u8NJuhpas8si9C3dsW1B4DHj0be6EnD8YoWy7nD0byWRjsepEgtnwyWNsCUfgMp69TlMO7R4q9NfuXJZw8A1b2od157eU7XvNpsVaRUbamog8xVUYfnoffZgmbomjoAAPcYNC+NSiEjnGtSnwC0JX/Dq7C1mDwGEADBDFCu8LbKbofKUwZQCOduPRwzKAnAcf0l6OTXV3IOxvUR58FB0SDdd/7W7qGYiGWWRjgqLfe627QaBPI0ED0siRZlojCvkT/qlbsCzcOF5qCRv7ZrDiBGXK7DXRfYtDn+KnuIseP/ZevW5tbi57aQfrjccdgEJT/pgVKzoLIsFxI1j9kgqmfn6OrepRCg6km5NiXmYGFAFPJjLB4pvBPWNQHHlCxPwtrybWLlUEgdjjQfHZjRQUfjfnVISKdRRVjJddvWDiLHXF4NAE0gIa/EukgPbILEt5Frkp/glqJGT3fUpeppQ3rfNZcGoNpm/2CsXW9C28GuoRxHbKztu9T+NoY1VoOD65GC5qDm/a5U9Oq51au2NbN5IC6fZvMQoUWmZB4h7qVO0s+mcfVXSmJqpzIXjcJahKVrSqDI78y2BIh73ouD+xpJD2yaxLehSKfee0YQ0xY576fyHriPdbxwu86M2tXbNt1XbhIakAmcC+m6gIVk5qAgZh4tPqtmsPRqmzQeoa2uZcJKO0iVkbJ+M82DYknW1ozsckEukn42jSvSK9hz72XAOTVmUZJI79vxqUZ69X3XjM2T+C5kWGgAgH9IOw1XSFedrn7hMljq19vCWhFoktoVoQO/pZ/0d6nX6VrZgKOM3VxN+vsS5Iru0UkAWHN+O7lFNFiOKYr8uSuAy/E79bfrboxhkCPLA/TjtPIKnAm7GisPvzwtoja9uSTVpUut6gPNpAE8U4BlnnPCNrjuW7FM9ThNlTaaIFj2FbEUifO4/TKfAWnK4oLPiKmKqAw/md/mYJBf/RXIryM9LpfrLCNc0HJ+bYMlgx4VUwIVxg0sdkX27c+gp03XlweV74W41KuwRv38wqv6njRzdTSVW8YF0uknoOm7y0jlJoi6dAHikrhtDIpMlbreYms3L2qaQwD3aoblXLSMdKLMP/iIGYWSoDj0oKosMMKk7Ta6oi4Bt7QXQ2e1U66tzZRVP0lDBAFXyUJUAI8req8x0GxNiThaEZ+IPktE3yGiN4noV4TrP09E3yyP3yOiHzh506xnrGpBtzukUemS9UlEr6VB9tgYpLKuektQXYgLA5h70JBUabtz6XN6/ZAb5GxqHfmC0LpHnj4g6QZH3wBAaWmPsQ8AmETFAOBya9W+bj4LawPAqqRXSGbFLklJGtYO8bs0dGaR/GsivEKjTCCiEMCvAfgpALcBfIOIvsbMr2vFvgvgv2Pmh0T0OQBfBvAj62tm2Rblp0MwO7l3rgx3Z3VK9SXqrxrUgLaDwSog9tYj+uSHggOJVeVSEN6zVIfP/uBdZkSxVKkv84m/qyJ/Q+h0flQY+eygLhUgVUt6mbffL7EVVPWMugcpU2l/oYXdSl8pWQPaSPwfBvAmM7/FzHMAXwXweaOdzL/HzA/Lj18H8Nz6mtgCHolXXG95zll2iZWCk9xjoUYcX3XV0qRwyUU2R33VANsk/RXY/Ns0PZMs+02kr9qmNAS77fYgNg+AaWj+nmp64PALoIyAmZwy3Ud625fAhu3Mo+/tCEIRqt3Uh8uyZ5mI41kA72qfb8Mvzf8mgH8tXSCiLwD4AgC88MILLZtowsulMjxXvhGLxJTOyv3PNjSOtjgJ+X3aigRy3yBJ/eIetJIkhpFQk1a1x7Pf7lnVoykpQULITiJNK38J4btPQn+2pjJbs51vwayf5GjJXuYkfc6EPCPDkccgvEK5pFyL5NMHn2UMqy3RRuK7zDn1gkQ/gYL4vyxdZ+YvM/PLzPzyjRv+PciWRhvdVAp2OaVH2eWb7lFlqkNPyuCzUywBQ3IbUxmTkFWRJilcSiD9+/nMFrbmoO6zn9O4+4/yr2ia96prKXnDWr3ekL7pkK1kzOruiWkWIJlGbtLr9TUNfmvoAwptJP5tAM9rn58D8L5diIi+H8BvAvgcM9+3r58aVlpSa1lvg9GtjTDXm+e7Z9XVwRoaMuz6JH8t1IHWo1lS5nbZrZb1lOdeC1HUKs2a9B1Vrkbp/UgZjn0edwxkB3EtvdrU8e45IwR67gXXkm2VXs02BLubsgraSPxvAPgEEb1ERD0APwvga3oBInoBwG8D+B+Z+Y31NrEBtlvZKqq1VEeLenxEdUl51znjfMsf2fn8pqY7vlse1b9+ExGZ3JpB0yuUlhspL99Hg9T3S2nHecbqpJdCtKVHCPYDpcKvLQHnGtAo8Zk5JaIvAfgdFImcv8LM3yaiL5bXXwHw9wBcA/DrVOg/KTO/fHrNXiM8zjJtJborONB3z6nB9vteg++C7QQE1N8Nh3ASTpL29l5+FfR3aadBs4uq1Fb6tabElZJzkXpnnp2ZavX4iiQByOFglKeyr4GxuWqb5DAnRCsXD2Z+FcCr1rlXtP//FoC/td6mrQG1DBiCLiugigFqacw7VSKvCt1oV1sdaDmqWeCgIDi50nJpKwCqjCK9rhVUpNfb1fSOlXpvl9O1A10o64NB299HHwDargRZCVWrog7yc+keTAy30XHNS3cSNtNl1wcfW1e0sJ8psU/yLNtYpOrSrMLO76It4SnJXrMFKPVfzcutAcBIahK1X6ID0Oi6W5TR3HebjKVKGxDVeHjyIqCMABWuSaGyNjRpzUkAGqUV2WsQ9ocQffa7IJ0WwmoTyO3COtqghF+TW7r+Im1FyFbtXQZJzS04KrNcNZHXG5STtahjmXfE5XKd2H6rolbRkijIqtenDwTS/P44qmV6MqCRv0q8oVyjRa1jPQv5G0d8CSutrwv3bzSsL18LvpHWgom9lus8lp1takY54iLrjlaXlOCzjcUeMKW/rXUU6+5lfS7pK/gX1PqGtV+dWAdZnyVk1OgrABTfQ1oyBQCkZBozpQF3nXEh2FDiL0NU36BwoQi/7rYwakk0KQCQme+iWot2zLd1knNgpvKqQUsqoebxagBoS/riZioI78gfUDXTNv7VbBnavz5N0beWau3Y5ESLaYBOfuM+FxirZZ1qgY0k/iowHFPWtWa+CbBURrHzk2Y0s2MUrBuynp351fMiA4aKQrWX39SAULPut+joEoGUX7s4p9ea2kh+KXOTalNTp/FkWq7akHh8CWptamr06niiie/6nS4c6dcyt3fP2dVlAEJSjTZ6sZ9QBpTUd5W3fBYCLXuSPQi02afQhnLBXQlNBKuuC6sLOnwZmgHTi1Aq2mS4XEN/uTgeBWvGhSP3WYJQU3MNuOaaRh3uixWx9L3/9ICeVYknNcOe/zoLav+6rPLr7BOsHeuoR//swbr69RMl8TeS7Cu2ue6YIxkxFpck8qtc8GLIbkMD2+xMbJd3ErhNVVJSEO8DIVrGxXn2SWEbAl1lmq55PMHW3befWIm/EVgX6X1okzU3QLOnsn2ehHN6ffYc3CY9W3/Jv9WX/exW6cJcU70W70QPmKrXK3xxX50+3wlX/dozTkOgdcR/EuBia3NMyAK+6+T429Ssplh+oUMH6eIwmiAY/dyJQpuZ4iK/RPbWxJPSqPtI34Q15HRwoSP+k4B1iYQ2Ur3pswQmwylHJJ1FkioDcHneSWa2jiVgk7wp6Mp8bsNgWysPs42nw+fW6Ij/pGJVN88GqS4KoXPqxKI2wfWlQxGW8bPN2Lk0+V04xfX5tuiI/ySgpTrYWmv0kF7/3x4ExLm9lcq77bOqy7lgSGQH6fX7zsPQa4R1N5Sx/2+odt14oqz6HdYAoZO1Wt62BgDAUulLJnKwcLkNlCqv7gsdz3P5rRuNgHcZc13kqbk5+Cq2rf2nmXB1SXQSfwPR2E88y0vGvS1U9hOlDAzg1aOVww6HDm3Abpt+2HDl/EfDd1hBM2itTTAVKn1+eka6VdER/7xwQlW0VT86oYA56dSg1TOaovl8/gJqEPD0YreLLupLikugkfxrILp3SfGE6Ih/Hlj3D+nrHW2W8QQpui7SNxE797jkukhfze897/EUV8JOFadJdh3dHP+scZo/qORn75P6jvXqVk48bbzVVNFwUWcwV+c0d19r3zhi94BQM+pJ7Vg1AesJBopl7AmusmdpkOyI/6RBkX+FTuyUpFJdLQ1wNgnyHjsdb/KeUIFKiU7NlvyiHVz9OW2J7wsCaxsGfl5u5h3xNxxiWsGgZa/Xyeoz+q0KIfmHmPPCNx0IHGVUvevY1uoUBgib/BctjqQj/llDdYY1dAQvtyW9XZLK6+r0rhUCaW/3JXeGESX9uvaxO0Wt4KKRXUdn3HvSoZJs6B55+tHk1y4Z/1Yhi93TWvQ8wyFIb+ZpuShfInQS/7zQxillbc9qUP21lFlFeVc9jv/1cz4jWWB9ZYeK7824e5HF6AahI/4TBmdIt2jx1wrpkn+dVrGWexkYRUTX3vaEb938CxIwcx7oiP+EonF3n5NKTp/fQNsElWvcseREY9UJl/I2ER3xn3CYLrqs/2mHJQcIJlSGt1VTey377LUoKJeM/B3xO9SxojZQ20+vSevQn3cR1r4uEfk74p8XztFGte704m2T01bPFB7OAa/mcHdJiLpudMS/pNDJuOpAsCzpivBdj29+Ww1hxed3WKAj/iWHnhNCoYl4p0043wDQkX096IjfoQZbAzhvsp3Z8y/RoNLKc4+IPktE3yGiN4noV4TrRET/pLz+TSL6ofU39QnDBe9kjem2z7ANHdaPRuITUQjg1wB8DsD3AfgbRPR9VrHPAfhEeXwBwG+suZ0dOpwuLtkA00bi/zCAN5n5LWaeA/gqgM9bZT4P4J9yga8D2Ceip9fc1icPl6yzdbg4aDPHfxbAu9rn2wB+pEWZZwF8oBcioi+g0AjwwgsvtGrgd7/0d1uV69ChQ3u0Sm0gnGsT3FmzyTLzl5n5ZWZ++caNG23a16FDh1NAG+LfBvC89vk5AO+vUKZDhw4XBG2I/w0AnyCil4ioB+BnAXzNKvM1AL9QWvd/FMBjZv7ArqhDhw4XA41zfGZOiehLAH4HRQT1V5j520T0xfL6KwBeBfDTAN4EMAbwi6fX5A4dOpwUrRx4mPlVFOTWz72i/c8Afmm9TevQocNpoUu91aHDJURH/A4dLiE64nfocAnREb9Dh0sIKuxy5/BgorsA3m5R9DqAe6fcnJPgIrfvIrcNuNjtu8htA9q370VmrnnLnRvx24KIXmPml8+7HS5c5PZd5LYBF7t9F7ltwMnb16n6HTpcQnTE79DhEmITiP/l825AAy5y+y5y24CL3b6L3DbghO278HP8Dh06rB+bIPE7dOiwZnTE79DhEuLCEP8iJ/Rs0bafL9v0TSL6PSL6gbNqW5v2aeU+Q0QZEf3MRWobEf04Ef0hEX2biP7dWbWtTfuIaI+I/iUR/VHZvjOLPCWirxDRHSL6luP66pxg5nM/UIT7/hmAjwPoAfgjAN9nlflpAP8aRbafHwXw+xeobX8ZwJXy/8+dVdvatk8r9/+iiLL8mYvSNgD7AF4H8EL5+eZFencA/mcA/7D8/waABwB6Z9S+/xbADwH4luP6ypy4KBL/Iif0bGwbM/8eMz8sP34dRQais0KbdwcAfxvAPwdw54K17ecA/DYzvwMAzHzR2scAdoiIAGyjIH56Fo1j5t8tn+fCypy4KMR3JetctsxpYNnn/k0Uo/BZobF9RPQsgL8O4BWcLdq8u08CuEJE/5aI/oCIfuHMWteufb8K4NMoUsn9MYC/w8z52TSvEStz4qLspLO2hJ6ngNbPJaKfQEH8HzvVFlmPFc7Z7fvHAH6ZmbNCcJ0Z2rQtAvCXAPwkgCGA/0hEX2fmN067cWjXvr8K4A8B/BUA3wPg/yaif8/MB6fctjZYmRMXhfgXOaFnq+cS0fcD+E0An2Pm+2fQLoU27XsZwFdL0l8H8NNElDLzv7gAbbsN4B4zHwM4JqLfBfADAM6C+G3a94sA/lcuJtVvEtF3AfwFAP/pDNrXhNU5cVaGlAYjRgTgLQAvYWFk+a+sMn8NpiHjP12gtr2AIt/gX76I784q/1s4O+Nem3f3aQD/T1l2BOBbAP7rC9S+3wDwv5T/3wLwHoDrZ/j7fgxu497KnLgQEp8vcELPlm37ewCuAfj1UqqmfEaRXS3bdy5o0zZm/hMi+jcAvgkgB/CbzCwuX51H+wD8AwC/RUR/jIJgv8zMZxKuS0T/DMCPA7hORLcB/H0Asda2lTnRuex26HAJcVGs+h06dDhDdMTv0OESoiN+hw6XEB3xO3S4hOiI36HDJURH/A4dLiE64nfocAnx/wNyldyN2jq+mAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 2\n",
    "print(\"Mode number {} of Karhunen-Loève decomposition\".format(n))\n",
    "show_mode(E[n, :], mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coordinates of degree of freedom number 300 are [0.27054928 0.42692123]\n"
     ]
    }
   ],
   "source": [
    "# cov = np.load(\"data/covariance.npy\", allow_pickle=True)\n",
    "# E = np.load(\"data/cov_modes.npy\", allow_pickle=True)\n",
    "\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "dofs = V.dofmap().dofs()\n",
    "\n",
    "# Get coordinates as len(dofs) x gdim array\n",
    "dim =V.dim()\n",
    "N = mesh.geometry().dim()\n",
    "dofs_x = V.tabulate_dof_coordinates()\n",
    "n_dof = 300\n",
    "print(\"Coordinates of degree of freedom number {0} are {1}\".format(n_dof, dofs_x[n_dof]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1668,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAEICAYAAAC9P1pMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARZklEQVR4nO3ce4xc5X3G8e+zO95dm2XBrAHD2o4NMhcTQQKO4+bSEpLUNm3q0KIKiEqDolpIcRqp/QPSPwItqIWokWgUwLKQ46YXO2pA1IlILKsNF4mbTQqGJTFZTMBr44BtYuPL2nv59Y85DsO8s96z67l48PORRtpz3nfOeTzsPHvOmTMoIjAzK9XS6ABmduJxMZhZwsVgZgkXg5klXAxmlnAxmFnCxfABJGm1pDtzzv2kpF9J2i/pizWOZk3CxWD/AHw3Ijoj4uHyQUnLJW2SdFjS6grjUyTdJ2mXpL2SHi8Zk6S7Je3OHt+SpJLx2ZJ+JumgpF9K+lyN/o02ToVGB7CG+xDQe4zxHcCdwCJgcoXxlRR/jy4G9gAfKRlbBnwRuAwIYAOwFViRja8BngKuzh4/lDQ3It6e2D/FqsVHDB8Akj4q6eeS3pX0A6CjbPyvJPVJ2iNpnaRzs/WvAucBP8pOJdrLtx0RD2VHErsr7PdC4E+AZRHxdkQMR8RzJVP+Evh2RPRHxHbg28CXs+deAFwO3BYRhyLiQeBF4M+O9/Ww4+diaHKS2oCHgX8DzgD+i5I3l6SrgH8C/hw4B3gdWAsQEecDbwBfyE4lDo9z9x/Ptvf32anEi5JK39iXAC+ULL+QrTs6tjUi3h1l3BrIxdD8FgKTgHsiYjAifghsLBn/ErAqIn6evfG/AfyepNlV2PcM4MPAXuBcYDnwr5IuzsY7s7Gj9gKd2XWG8rGj46dWIZcdJxdD8zsX2B7v/zbc62Xjv1uOiP0UTwt6qrDvQ8AgcGdEHImIx4CfAX+Yje8HukrmdwH7s6zlY0fH38UazsXQ/N4Eekqv9gOzSn7eQfECIwCSTgG6ge1V2PfmMcZ7KV54POoy3rvQ2QucJ+nUUcatgVwMze8pYAj4a0kFSX8KLCgZ/0/gJkkfyS4u/iPwTET8Os/Gs212AK1Aq6QOSUc/zXqc4jWKb2TzPglcCazPxr8P/I2knuyC598CqwEi4hXgeeC2bJvXAJcCD07kRbAqiwg/mvwBzAf+j+Jh+A+yx50l4zcDr1L8OPHHwIySsV8DnzvGtm+n+FFj6eP2kvFLKJbTAeBl4JqSMQHfyva7J/tZJeOzgUcpnpJsOVYOP+r7UPYfyMzsd3wqYWaJMYtB0ipJb0l6aZRxSfpOdgPNZkmXVz+mmdVTniOG1cDiY4wvAeZmj2XA/ccfy8waacxiiIjHKV44Gs1S4PtR9DRwuqRzqhXQzOqvGl+i6gG2lSz3Z+veLJ8oaRnFowpaab1iik6rwu7NbDTvxp5dEXHmeJ9XjWJQhXUVP+qIiJUUv41HV0t3LCwsqsLuzWw0GwbXvD72rFQ1PpXoB2aWLM+geLedmTWpahTDOuDG7NOJhcDeiEhOI8yseYx5KiFpDcXbXKdJ6gduo/htPiJiBfAIxf/JRh9wELipVmHNrD7GLIaIuH6M8QC+WrVEZtZwvvPRzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgazE1Brz3Rau6c2bP+Fhu3ZzCrSJXPZO7eL4XbRejgoHB5h0r4hWgeG0MAQLQNHYOAww9t31ixDrmKQtBj4F6AVeCAi7iobPw34d2BWts1/jojvVTmr2Qdea890DszqZN/sFg7MGgGgsL9Ay2CB1sPQcgQ69gRn9B6ARhaDpFbgXuDzQD+wUdK6iHi5ZNpXgZcj4guSzgS2SPqPiDhSk9RmH0Ct3VM5dOHZvNtT4ND0YNL0g8w9620ADgy2ceBIG29vmwp9BQbO6uCUnuk1O2rIc8SwAOiLiK0AktYCS4HSYgjgVEkCOoE9wFCVs5p9YLVM7mB49nT2z5jE/lkQ0wf4aM92/nja8yzoeIPdIx2sfvvTPHpgMofP6OTgma10TJ9as6OGPBcfe4BtJcv92bpS3wUuBnYALwJfj4iR8g1JWiZpk6RNgzEwwchmH0zDHQWG2wRAoW2Y0yYdoruwn1mFds4rDPChybs5d+peBs8YYuAMMXBWB60902uSJc8Rgyqsi7LlRcDzwFXA+cAGSU9ExL73PSliJbASoKulu3wbZietkUMDtP2yn9PaZzHS1sa+zik81TabrkkDQC+7h87m9UPdHDjSBu0jDJ4Kh7uKRw0te37LyKHq/qHNUwz9wMyS5RkUjwxK3QTcFREB9El6DbgIeLYqKc1OAsO732HylnZG2s8BCuyji/VcxLYzix9bvvz22by7s5P23xSY8ia07xuuWZY8xbARmCtpDrAduA64oWzOG8BngScknQ1cCGytZlCzk8Hw9p2cAkCxHA7uP51nzplCy6RhYk87U3a2MGVn0Nk/SPtvDqKt26p+tAA5iiEihiQtB9ZT/LhyVUT0Sro5G18B3AGslvQixVOPWyJiV9XTmp0Ehrfv5JSOdloOn86k/ZM4uL+NwU5o3wOnvBl09h+mbec+hl99PTmnrxYVj/7rr6ulOxYWFjVk32bNoGVyB3HeTA7NKn4K0bYv6Nx2kNad7+T+mHLD4JrnImL+ePftOx/NTlAjhwag91ec8tvpdEyfWrzrcecuhne/U/N9uxjMTnDD23fC9p0EULvLje/nL1GZWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlnAxmFnCxWBmCReDmSVcDGaWcDGYWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlnAxmFnCxWBmCReDmSVcDGaWcDGYWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlshVDJIWS9oiqU/SraPMuVLS85J6JT1W3ZhmVk+FsSZIagXuBT4P9AMbJa2LiJdL5pwO3Acsjog3JJ1Vo7xmVgd5jhgWAH0RsTUijgBrgaVlc24AHoqINwAi4q3qxjSzespTDD3AtpLl/mxdqQuAqZIelfScpBsrbUjSMkmbJG0ajIGJJTazmhvzVAJQhXVRYTtXAJ8FJgNPSXo6Il5535MiVgIrAbpausu3YWYniDzF0A/MLFmeAeyoMGdXRBwADkh6HLgMeAUzazp5TiU2AnMlzZHUBlwHrCub89/ApyUVJE0BPg78orpRzaxexjxiiIghScuB9UArsCoieiXdnI2viIhfSPopsBkYAR6IiJdqGdzMakcRjTnV72rpjoWFRQ3Zt9nJYsPgmuciYv54n+c7H80s4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMErmKQdJiSVsk9Um69RjzPiZpWNK11YtoZvU2ZjFIagXuBZYA84DrJc0bZd7dwPpqhzSz+spzxLAA6IuIrRFxBFgLLK0w72vAg8BbVcxnZg2Qpxh6gG0ly/3Zut+R1ANcA6w41oYkLZO0SdKmwRgYb1Yzq5M8xaAK66Js+R7glogYPtaGImJlRMyPiPmT1JEzopnVWyHHnH5gZsnyDGBH2Zz5wFpJANOAqyUNRcTD1QhpZvWVpxg2AnMlzQG2A9cBN5ROiIg5R3+WtBr4sUvBrHmNWQwRMSRpOcVPG1qBVRHRK+nmbPyY1xXMrPnkOWIgIh4BHilbV7EQIuLLxx/LzBrJdz6aWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlnAxmFnCxWBmCReDmSVcDGaWcDGYWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlnAxmFnCxWBmCReDmSVcDGaWcDGYWcLFYGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlnAxmFnCxWBmiVzFIGmxpC2S+iTdWmH8S5I2Z48nJV1W/ahmVi9jFoOkVuBeYAkwD7he0ryyaa8BfxARlwJ3ACurHdTM6ifPEcMCoC8itkbEEWAtsLR0QkQ8GRHvZItPAzOqG9PM6ilPMfQA20qW+7N1o/kK8JNKA5KWSdokadNgDORPaWZ1VcgxRxXWRcWJ0mcoFsOnKo1HxEqy04yulu6K2zCzxstTDP3AzJLlGcCO8kmSLgUeAJZExO7qxDOzRshzKrERmCtpjqQ24DpgXekESbOAh4C/iIhXqh/TzOppzCOGiBiStBxYD7QCqyKiV9LN2fgK4JtAN3CfJIChiJhfu9hmVkuKaMypfldLdywsLGrIvs1OFhsG1zw3kT/SvvPRzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCzhYjCzhIvBzBIuBjNLuBjMLOFiMLOEi8HMEi4GM0u4GMws4WIws4SLwcwSLgYzS7gYzCyRqxgkLZa0RVKfpFsrjEvSd7LxzZIur35UM6uXMYtBUitwL7AEmAdcL2le2bQlwNzssQy4v8o5zayO8hwxLAD6ImJrRBwB1gJLy+YsBb4fRU8Dp0s6p8pZzaxOCjnm9ADbSpb7gY/nmNMDvFk6SdIyikcUAIc3DK55aVxpG2sasKvRIXJqpqzQXHmbKSvAhRN5Up5iUIV1MYE5RMRKYCWApE0RMT/H/k8IzZS3mbJCc+VtpqxQzDuR5+U5legHZpYszwB2TGCOmTWJPMWwEZgraY6kNuA6YF3ZnHXAjdmnEwuBvRHxZvmGzKw5jHkqERFDkpYD64FWYFVE9Eq6ORtfATwCXA30AQeBm3Lse+WEUzdGM+VtpqzQXHmbKStMMK8ikksBZnaS852PZpZwMZhZoubF0Ey3U+fI+qUs42ZJT0q6rBE5S/IcM2/JvI9JGpZ0bT3zlWUYM6ukKyU9L6lX0mP1zliWZazfhdMk/UjSC1nePNfVakLSKklvSap4X9CE3mMRUbMHxYuVrwLnAW3AC8C8sjlXAz+heC/EQuCZWmY6zqyfAKZmPy9pVNa8eUvm/S/FC8TXnqhZgdOBl4FZ2fJZJ/JrC/wdcHf285nAHqCtQXl/H7gceGmU8XG/x2p9xNBMt1OPmTUinoyId7LFpyner9EoeV5bgK8BDwJv1TNcmTxZbwAeiog3ACLiRM8bwKmSBHRSLIah+sbMgkQ8nu1/NON+j9W6GEa7VXq8c+phvDm+QrGFG2XMvJJ6gGuAFXXMVUme1/YCYKqkRyU9J+nGuqVL5cn7XeBiijfyvQh8PSJG6hNv3Mb9HstzS/TxqNrt1HWQO4ekz1Ashk/VNNGx5cl7D3BLRAwX/7A1TJ6sBeAK4LPAZOApSU9HxCu1DldBnryLgOeBq4DzgQ2SnoiIfTXONhHjfo/Vuhia6XbqXDkkXQo8ACyJiN11ylZJnrzzgbVZKUwDrpY0FBEP1yXhe/L+HuyKiAPAAUmPA5cBjSiGPHlvAu6K4kl8n6TXgIuAZ+sTcVzG/x6r8UWRArAVmMN7F3EuKZvzR7z/wsizDbqAkyfrLIp3d36iERnHm7ds/moad/Exz2t7MfA/2dwpwEvAh0/gvPcDt2c/nw1sB6Y18PdhNqNffBz3e6ymRwxRu9upG5X1m0A3cF/2V3goGvRNu5x5Twh5skbELyT9FNgMjAAPRERDvpaf87W9A1gt6UWKb7hbIqIhX8eWtAa4EpgmqR+4DZhUknXc7zHfEm1mCd/5aGYJF4OZJVwMZpZwMZhZwsVgZgkXg5klXAxmlvh/44qs2VVRKPQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh = Mesh(\"data/mesh.xml\")\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "u = Function(V)\n",
    "print(np.array(u.vector()[:]).shape)\n",
    "\n",
    "for j in range(16):\n",
    "    for i in range(1668):\n",
    "        if i == (j+1)*100:\n",
    "            u.vector()[i]= 1\n",
    "        else:\n",
    "            u.vector()[i]=0\n",
    "    plot(u, title='dof {}'.format((1+j)*100))\n",
    "    plt.savefig('data/component_{}.png'.format((1+j)*100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution number : 0\n",
      "Solution number : 1\n",
      "Solution number : 2\n",
      "Solution number : 3\n",
      "Solution number : 4\n",
      "Solution number : 5\n",
      "Solution number : 6\n",
      "Solution number : 7\n",
      "Solution number : 8\n",
      "Solution number : 9\n",
      "Solution number : 10\n",
      "Solution number : 11\n",
      "Solution number : 12\n",
      "Solution number : 13\n",
      "Solution number : 14\n",
      "Solution number : 15\n",
      "Solution number : 16\n",
      "Solution number : 17\n",
      "Solution number : 18\n",
      "Solution number : 19\n",
      "Solution number : 20\n",
      "Solution number : 21\n",
      "Solution number : 22\n",
      "Solution number : 23\n",
      "Solution number : 24\n",
      "Solution number : 25\n",
      "Solution number : 26\n",
      "Solution number : 27\n",
      "Solution number : 28\n",
      "Solution number : 29\n",
      "Solution number : 30\n",
      "Solution number : 31\n",
      "Solution number : 32\n",
      "Solution number : 33\n",
      "Solution number : 34\n",
      "Solution number : 35\n",
      "Solution number : 36\n",
      "Solution number : 37\n",
      "Solution number : 38\n",
      "Solution number : 39\n",
      "Solution number : 40\n",
      "Solution number : 41\n",
      "Solution number : 42\n",
      "Solution number : 43\n",
      "Solution number : 44\n",
      "Solution number : 45\n",
      "Solution number : 46\n",
      "Solution number : 47\n",
      "Solution number : 48\n",
      "Solution number : 49\n",
      "Solution number : 50\n",
      "Solution number : 51\n",
      "Solution number : 52\n",
      "Solution number : 53\n",
      "Solution number : 54\n",
      "Solution number : 55\n",
      "Solution number : 56\n",
      "Solution number : 57\n",
      "Solution number : 58\n",
      "Solution number : 59\n",
      "Solution number : 60\n",
      "Solution number : 61\n",
      "Solution number : 62\n",
      "Solution number : 63\n",
      "Solution number : 64\n",
      "Solution number : 65\n",
      "Solution number : 66\n",
      "Solution number : 67\n",
      "Solution number : 68\n",
      "Solution number : 69\n",
      "Solution number : 70\n",
      "Solution number : 71\n",
      "Solution number : 72\n",
      "Solution number : 73\n",
      "Solution number : 74\n",
      "Solution number : 75\n",
      "Solution number : 76\n",
      "Solution number : 77\n",
      "Solution number : 78\n",
      "Solution number : 79\n",
      "Solution number : 80\n",
      "Solution number : 81\n",
      "Solution number : 82\n",
      "Solution number : 83\n",
      "Solution number : 84\n",
      "Solution number : 85\n",
      "Solution number : 86\n",
      "Solution number : 87\n",
      "Solution number : 88\n",
      "Solution number : 89\n",
      "Solution number : 90\n",
      "Solution number : 91\n",
      "Solution number : 92\n",
      "Solution number : 93\n",
      "Solution number : 94\n",
      "Solution number : 95\n",
      "Solution number : 96\n",
      "Solution number : 97\n",
      "Solution number : 98\n",
      "Solution number : 99\n",
      "Solution number : 100\n",
      "Solution number : 101\n",
      "Solution number : 102\n",
      "Solution number : 103\n",
      "Solution number : 104\n",
      "Solution number : 105\n",
      "Solution number : 106\n",
      "Solution number : 107\n",
      "Solution number : 108\n",
      "Solution number : 109\n",
      "Solution number : 110\n",
      "Solution number : 111\n",
      "Solution number : 112\n",
      "Solution number : 113\n",
      "Solution number : 114\n",
      "Solution number : 115\n",
      "Solution number : 116\n",
      "Solution number : 117\n",
      "Solution number : 118\n",
      "Solution number : 119\n",
      "Solution number : 120\n",
      "Solution number : 121\n",
      "Solution number : 122\n",
      "Solution number : 123\n",
      "Solution number : 124\n",
      "Solution number : 125\n",
      "Solution number : 126\n",
      "Solution number : 127\n",
      "Solution number : 128\n",
      "Solution number : 129\n",
      "Solution number : 130\n",
      "Solution number : 131\n",
      "Solution number : 132\n",
      "Solution number : 133\n",
      "Solution number : 134\n",
      "Solution number : 135\n",
      "Solution number : 136\n",
      "Solution number : 137\n",
      "Solution number : 138\n",
      "Solution number : 139\n",
      "Solution number : 140\n",
      "Solution number : 141\n",
      "Solution number : 142\n",
      "Solution number : 143\n",
      "Solution number : 144\n",
      "Solution number : 145\n",
      "Solution number : 146\n",
      "Solution number : 147\n",
      "Solution number : 148\n",
      "Solution number : 149\n",
      "Solution number : 150\n",
      "Solution number : 151\n",
      "Solution number : 152\n",
      "Solution number : 153\n",
      "Solution number : 154\n",
      "Solution number : 155\n",
      "Solution number : 156\n",
      "Solution number : 157\n",
      "Solution number : 158\n",
      "Solution number : 159\n",
      "Solution number : 160\n",
      "Solution number : 161\n",
      "Solution number : 162\n",
      "Solution number : 163\n",
      "Solution number : 164\n",
      "Solution number : 165\n",
      "Solution number : 166\n",
      "Solution number : 167\n",
      "Solution number : 168\n",
      "Solution number : 169\n",
      "Solution number : 170\n",
      "Solution number : 171\n",
      "Solution number : 172\n",
      "Solution number : 173\n",
      "Solution number : 174\n",
      "Solution number : 175\n",
      "Solution number : 176\n",
      "Solution number : 177\n",
      "Solution number : 178\n",
      "Solution number : 179\n",
      "Solution number : 180\n",
      "Solution number : 181\n",
      "Solution number : 182\n",
      "Solution number : 183\n",
      "Solution number : 184\n",
      "Solution number : 185\n",
      "Solution number : 186\n",
      "Solution number : 187\n",
      "Solution number : 188\n",
      "Solution number : 189\n",
      "Solution number : 190\n",
      "Solution number : 191\n",
      "Solution number : 192\n",
      "Solution number : 193\n",
      "Solution number : 194\n",
      "Solution number : 195\n",
      "Solution number : 196\n",
      "Solution number : 197\n",
      "Solution number : 198\n",
      "Solution number : 199\n",
      "Solution number : 200\n",
      "Solution number : 201\n",
      "Solution number : 202\n",
      "Solution number : 203\n",
      "Solution number : 204\n",
      "Solution number : 205\n",
      "Solution number : 206\n",
      "Solution number : 207\n",
      "Solution number : 208\n",
      "Solution number : 209\n",
      "Solution number : 210\n",
      "Solution number : 211\n",
      "Solution number : 212\n",
      "Solution number : 213\n",
      "Solution number : 214\n",
      "Solution number : 215\n",
      "Solution number : 216\n",
      "Solution number : 217\n",
      "Solution number : 218\n",
      "Solution number : 219\n",
      "Solution number : 220\n",
      "Solution number : 221\n",
      "Solution number : 222\n",
      "Solution number : 223\n",
      "Solution number : 224\n",
      "Solution number : 225\n",
      "Solution number : 226\n",
      "Solution number : 227\n",
      "Solution number : 228\n",
      "Solution number : 229\n",
      "Solution number : 230\n",
      "Solution number : 231\n",
      "Solution number : 232\n",
      "Solution number : 233\n",
      "Solution number : 234\n",
      "Solution number : 235\n",
      "Solution number : 236\n",
      "Solution number : 237\n",
      "Solution number : 238\n",
      "Solution number : 239\n",
      "Solution number : 240\n",
      "Solution number : 241\n",
      "Solution number : 242\n",
      "Solution number : 243\n",
      "Solution number : 244\n",
      "Solution number : 245\n",
      "Solution number : 246\n",
      "Solution number : 247\n",
      "Solution number : 248\n",
      "Solution number : 249\n",
      "Solution number : 250\n",
      "Solution number : 251\n",
      "Solution number : 252\n",
      "Solution number : 253\n",
      "Solution number : 254\n",
      "Solution number : 255\n",
      "Solution number : 256\n",
      "Solution number : 257\n",
      "Solution number : 258\n",
      "Solution number : 259\n",
      "Solution number : 260\n",
      "Solution number : 261\n",
      "Solution number : 262\n",
      "Solution number : 263\n",
      "Solution number : 264\n",
      "Solution number : 265\n",
      "Solution number : 266\n",
      "Solution number : 267\n",
      "Solution number : 268\n",
      "Solution number : 269\n",
      "Solution number : 270\n",
      "Solution number : 271\n",
      "Solution number : 272\n",
      "Solution number : 273\n",
      "Solution number : 274\n",
      "Solution number : 275\n",
      "Solution number : 276\n",
      "Solution number : 277\n",
      "Solution number : 278\n",
      "Solution number : 279\n",
      "Solution number : 280\n",
      "Solution number : 281\n",
      "Solution number : 282\n",
      "Solution number : 283\n",
      "Solution number : 284\n",
      "Solution number : 285\n",
      "Solution number : 286\n",
      "Solution number : 287\n",
      "Solution number : 288\n",
      "Solution number : 289\n",
      "Solution number : 290\n",
      "Solution number : 291\n",
      "Solution number : 292\n",
      "Solution number : 293\n",
      "Solution number : 294\n",
      "Solution number : 295\n",
      "Solution number : 296\n",
      "Solution number : 297\n",
      "Solution number : 298\n",
      "Solution number : 299\n",
      "Solution number : 300\n",
      "Solution number : 301\n",
      "Solution number : 302\n",
      "Solution number : 303\n",
      "Solution number : 304\n",
      "Solution number : 305\n",
      "Solution number : 306\n",
      "Solution number : 307\n",
      "Solution number : 308\n",
      "Solution number : 309\n",
      "Solution number : 310\n",
      "Solution number : 311\n",
      "Solution number : 312\n",
      "Solution number : 313\n",
      "Solution number : 314\n",
      "Solution number : 315\n",
      "Solution number : 316\n",
      "Solution number : 317\n",
      "Solution number : 318\n",
      "Solution number : 319\n",
      "Solution number : 320\n",
      "Solution number : 321\n",
      "Solution number : 322\n",
      "Solution number : 323\n",
      "Solution number : 324\n",
      "Solution number : 325\n",
      "Solution number : 326\n",
      "Solution number : 327\n",
      "Solution number : 328\n",
      "Solution number : 329\n",
      "Solution number : 330\n",
      "Solution number : 331\n",
      "Solution number : 332\n",
      "Solution number : 333\n",
      "Solution number : 334\n",
      "Solution number : 335\n",
      "Solution number : 336\n",
      "Solution number : 337\n",
      "Solution number : 338\n",
      "Solution number : 339\n",
      "Solution number : 340\n",
      "Solution number : 341\n",
      "Solution number : 342\n",
      "Solution number : 343\n",
      "Solution number : 344\n",
      "Solution number : 345\n",
      "Solution number : 346\n",
      "Solution number : 347\n",
      "Solution number : 348\n",
      "Solution number : 349\n",
      "Solution number : 350\n",
      "Solution number : 351\n",
      "Solution number : 352\n",
      "Solution number : 353\n",
      "Solution number : 354\n",
      "Solution number : 355\n",
      "Solution number : 356\n",
      "Solution number : 357\n",
      "Solution number : 358\n",
      "Solution number : 359\n",
      "Solution number : 360\n",
      "Solution number : 361\n",
      "Solution number : 362\n",
      "Solution number : 363\n",
      "Solution number : 364\n",
      "Solution number : 365\n",
      "Solution number : 366\n",
      "Solution number : 367\n",
      "Solution number : 368\n",
      "Solution number : 369\n",
      "Solution number : 370\n",
      "Solution number : 371\n",
      "Solution number : 372\n",
      "Solution number : 373\n",
      "Solution number : 374\n",
      "Solution number : 375\n",
      "Solution number : 376\n",
      "Solution number : 377\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution number : 378\n",
      "Solution number : 379\n",
      "Solution number : 380\n",
      "Solution number : 381\n",
      "Solution number : 382\n",
      "Solution number : 383\n",
      "Solution number : 384\n",
      "Solution number : 385\n",
      "Solution number : 386\n",
      "Solution number : 387\n",
      "Solution number : 388\n",
      "Solution number : 389\n",
      "Solution number : 390\n",
      "Solution number : 391\n",
      "Solution number : 392\n",
      "Solution number : 393\n",
      "Solution number : 394\n",
      "Solution number : 395\n",
      "Solution number : 396\n",
      "Solution number : 397\n",
      "Solution number : 398\n",
      "Solution number : 399\n",
      "Solution number : 400\n",
      "Solution number : 401\n",
      "Solution number : 402\n",
      "Solution number : 403\n",
      "Solution number : 404\n",
      "Solution number : 405\n",
      "Solution number : 406\n",
      "Solution number : 407\n",
      "Solution number : 408\n",
      "Solution number : 409\n",
      "Solution number : 410\n",
      "Solution number : 411\n",
      "Solution number : 412\n",
      "Solution number : 413\n",
      "Solution number : 414\n",
      "Solution number : 415\n",
      "Solution number : 416\n",
      "Solution number : 417\n",
      "Solution number : 418\n",
      "Solution number : 419\n",
      "Solution number : 420\n",
      "Solution number : 421\n",
      "Solution number : 422\n",
      "Solution number : 423\n",
      "Solution number : 424\n",
      "Solution number : 425\n",
      "Solution number : 426\n",
      "Solution number : 427\n",
      "Solution number : 428\n",
      "Solution number : 429\n",
      "Solution number : 430\n",
      "Solution number : 431\n",
      "Solution number : 432\n",
      "Solution number : 433\n",
      "Solution number : 434\n",
      "Solution number : 435\n",
      "Solution number : 436\n",
      "Solution number : 437\n",
      "Solution number : 438\n",
      "Solution number : 439\n",
      "Solution number : 440\n",
      "Solution number : 441\n",
      "Solution number : 442\n",
      "Solution number : 443\n",
      "Solution number : 444\n",
      "Solution number : 445\n",
      "Solution number : 446\n",
      "Solution number : 447\n",
      "Solution number : 448\n",
      "Solution number : 449\n",
      "Solution number : 450\n",
      "Solution number : 451\n",
      "Solution number : 452\n",
      "Solution number : 453\n",
      "Solution number : 454\n",
      "Solution number : 455\n",
      "Solution number : 456\n",
      "Solution number : 457\n",
      "Solution number : 458\n",
      "Solution number : 459\n",
      "Solution number : 460\n",
      "Solution number : 461\n",
      "Solution number : 462\n",
      "Solution number : 463\n",
      "Solution number : 464\n",
      "Solution number : 465\n",
      "Solution number : 466\n",
      "Solution number : 467\n",
      "Solution number : 468\n",
      "Solution number : 469\n",
      "Solution number : 470\n",
      "Solution number : 471\n",
      "Solution number : 472\n",
      "Solution number : 473\n",
      "Solution number : 474\n",
      "Solution number : 475\n",
      "Solution number : 476\n",
      "Solution number : 477\n",
      "Solution number : 478\n",
      "Solution number : 479\n",
      "Solution number : 480\n",
      "Solution number : 481\n",
      "Solution number : 482\n",
      "Solution number : 483\n",
      "Solution number : 484\n",
      "Solution number : 485\n",
      "Solution number : 486\n",
      "Solution number : 487\n",
      "Solution number : 488\n",
      "Solution number : 489\n",
      "Solution number : 490\n",
      "Solution number : 491\n",
      "Solution number : 492\n",
      "Solution number : 493\n",
      "Solution number : 494\n",
      "Solution number : 495\n",
      "Solution number : 496\n",
      "Solution number : 497\n",
      "Solution number : 498\n",
      "Solution number : 499\n"
     ]
    }
   ],
   "source": [
    "for it in range(M):\n",
    "  print(\"Solution number :\", it)\n",
    "  #set conductivity\n",
    "  kappa = set_conductivity(it, mesh, c)\n",
    "  #plot(kappa)\n",
    "  #plt.show()\n",
    "\n",
    "  #compute solution\n",
    "  u = compute_solution(it, mesh, kappa, pl=False)#pl=True to plot\n",
    "  u_res = restrict(mesh, u)\n",
    "  #print(\"mean of the solution restricted on the outflow (right side)\", u_res)\n",
    "  \n",
    "  #compute gradients\n",
    "  for j in range(m):\n",
    "    #print(\"Evaluating gradient component number :\", j)\n",
    "    du = compute_gradients(j, mesh, kappa, E, boundary, cache, u)\n",
    "    du_res = restrict(mesh, du)\n",
    "    cache[:, j] = du.vector()[:]\n",
    "    cache_res[j] = du_res\n",
    "\n",
    "  file = Path(\"data/outputs.npy\")\n",
    "  with file.open('ab') as f:\n",
    "      np.save(f, u.vector()[:])\n",
    "  \n",
    "  file = Path(\"data/outputs_res.npy\")\n",
    "  with file.open('ab') as f:\n",
    "      np.save(f, u_res)\n",
    "\n",
    "  file = Path(\"data/gradients.npy\")\n",
    "  with file.open('ab') as f:\n",
    "      np.save(f, cache)\n",
    "\n",
    "  file = Path(\"data/gradients_res.npy\")\n",
    "  with file.open('ab') as f:\n",
    "      np.save(f, cache_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks\n",
    "\n",
    "Try different correlation coefficients beta=0.03, 0.04 and compare Karhunen-Loève modes with vectorial active subspaces modes in 05_SPDE_on_fenics_modes.ipynb\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "athena",
   "language": "python",
   "name": "athena"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}