laplace-fvm-fracture.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:29.867150Z",
     "iopub.status.busy": "2023-10-06T07:42:29.866386Z",
     "iopub.status.idle": "2023-10-06T07:42:30.284642Z",
     "shell.execute_reply": "2023-10-06T07:42:30.285164Z",
     "shell.execute_reply.started": "2023-10-06T07:34:33.107333Z"
    },
    "papermill": {
     "duration": 0.431143,
     "end_time": "2023-10-06T07:42:30.285397",
     "exception": false,
     "start_time": "2023-10-06T07:42:29.854254",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "from scipy.sparse import csc_matrix, linalg\n",
    "\n",
    "from math import factorial, floor\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "font = {'size'   : 14}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:30.305670Z",
     "iopub.status.busy": "2023-10-06T07:42:30.304988Z",
     "iopub.status.idle": "2023-10-06T07:42:30.314017Z",
     "shell.execute_reply": "2023-10-06T07:42:30.313224Z",
     "shell.execute_reply.started": "2023-10-06T07:34:33.464104Z"
    },
    "papermill": {
     "duration": 0.019076,
     "end_time": "2023-10-06T07:42:30.314156",
     "exception": false,
     "start_time": "2023-10-06T07:42:30.295080",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def Vfunction(i, N):\n",
    "    sum1 = 0.\n",
    "    kmin = int(floor((i + 1.) / 2.))\n",
    "    kmax = min(i, N)\n",
    "    for k in range(kmin, kmax+1):\n",
    "        sum1 += (1.*(k**N) * factorial(2 * k) / (factorial(N - k) * factorial(k) * factorial(k - 1) * factorial(i - k) * factorial(2 * k - i)))\n",
    "    return ((-1.)**(N + i)) * sum1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:30.334680Z",
     "iopub.status.busy": "2023-10-06T07:42:30.334031Z",
     "iopub.status.idle": "2023-10-06T07:42:30.336895Z",
     "shell.execute_reply": "2023-10-06T07:42:30.337347Z",
     "shell.execute_reply.started": "2023-10-06T07:34:33.473105Z"
    },
    "papermill": {
     "duration": 0.016218,
     "end_time": "2023-10-06T07:42:30.337531",
     "exception": false,
     "start_time": "2023-10-06T07:42:30.321313",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def solve(matrixR,vectorF):\n",
    "    return linalg.spsolve(csc_matrix(matrixR, dtype=float), csc_matrix(vectorF, dtype=float)) # only for nCells <= 100 x 100\n",
    "\n",
    "    # Iterative solvers###################################\n",
    "    # Minimal residual iterative method (krylov) # can reach to nCells = 400 x 400 or more\n",
    "    '''\n",
    "    A = csc_matrix(matrixR, dtype=float)\n",
    "    A = A + A.T\n",
    "    vectorU, exitCode = scipy.sparse.linalg.minres(A, vectorF,tol=1e-100, maxiter=10000, show=False)\n",
    "    print(np.allclose(A.dot(vectorF), vectorU))\n",
    "    return vectorU'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:30.355408Z",
     "iopub.status.busy": "2023-10-06T07:42:30.354761Z",
     "iopub.status.idle": "2023-10-06T07:42:30.396001Z",
     "shell.execute_reply": "2023-10-06T07:42:30.396537Z",
     "shell.execute_reply.started": "2023-10-06T07:34:33.487671Z"
    },
    "papermill": {
     "duration": 0.052037,
     "end_time": "2023-10-06T07:42:30.396708",
     "exception": false,
     "start_time": "2023-10-06T07:42:30.344671",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rigidity(n, s, x,y, cxx,cyy,cf, phi,e):\n",
    "    nCell = n*n\n",
    "    M = np.zeros((nCell,nCell))\n",
    "    for i in range(1,n-1):\n",
    "        for j in range(1,n-1):\n",
    "            iCell = i*n + j\n",
    "            iW = (i-1)*n + j\n",
    "            iE = (i+1)*n + j\n",
    "            iS = iCell - 1\n",
    "            iN = iCell + 1\n",
    "\n",
    "            l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "            l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "            c1xi  = (cxx[i][j] + cxx[i-1][j]) /2. /l2 /(x[i][j] - x[i-1][j])\n",
    "            c2eta =-(cyy[i][j+1] + cyy[i][j]) /2. /l1 /(y[i][j+1] - y[i][j])\n",
    "            c3xi  =-(cxx[i+1][j] + cxx[i][j]) /2. /l2 /(x[i+1][j] - x[i][j])\n",
    "            c4eta = (cyy[i][j] + cyy[i][j-1]) /2. /l1 /(y[i][j] - y[i][j-1])\n",
    "\n",
    "            M[iCell][iCell] = s + c1xi - c2eta - c3xi + c4eta\n",
    "\n",
    "            M[iCell][iW] = -c1xi         \n",
    "            M[iCell][iN] = c2eta\n",
    "            M[iCell][iE] = c3xi\n",
    "            M[iCell][iS] = -c4eta\n",
    "\n",
    "\n",
    "    # Fracture\n",
    "    i = int(n/2.)\n",
    "    for j in range(int(40.*n/101.),int(60.*n/101.)): \n",
    "    # the end range of int(3.*n/4.)-1 caused a day long issue\n",
    "\n",
    "        iCell = i*n + j\n",
    "        for k in range(0,nCell):\n",
    "            M[iCell][k] = 0. \n",
    "\n",
    "\n",
    "        iW = (i-1)*n + j\n",
    "        iE = (i+1)*n + j\n",
    "        iS = iCell - 1\n",
    "        iN = iCell + 1\n",
    "\n",
    "        l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "        l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "        cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "\n",
    "        cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "\n",
    "        M[iCell][iCell] = s + 2.*cf/l1**2. + cfW - cfE\n",
    "\n",
    "\n",
    "\n",
    "        M[iCell][iN] = -cf/l1**2.\n",
    "        M[iCell][iS] = -cf/l1**2.\n",
    "        M[iCell][iW] = -cfW         \n",
    "        M[iCell][iE] = cfE\n",
    "\n",
    "\n",
    "\n",
    "    # Fracture tip 1\n",
    "    i = int(n/2.)\n",
    "    j = int(40.*n/101.)\n",
    "\n",
    "\n",
    "    iCell = i*n + j\n",
    "    for k in range(0,nCell):\n",
    "        M[iCell][k] = 0. \n",
    "\n",
    "    iW = (i-1)*n + j\n",
    "    iE = (i+1)*n + j\n",
    "    iS = iCell - 1\n",
    "    iN = iCell + 1\n",
    "\n",
    "    l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "    l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "\n",
    "\n",
    "    cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "\n",
    "    cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "\n",
    "    # Removing ctS allow a best fit for high matrix permeability, steady state regime\n",
    "    #coeffff = 0.\n",
    "    #ctS = coeffff*phi*l2/e*(cyy[i][j] + cyy[i][j-1]) /2. /(y[i][j] - y[i][j-1])\n",
    "\n",
    "    M[iCell][iCell] = s + cf/l1**2. + cfW - cfE #+ ctS\n",
    "\n",
    "\n",
    "    M[iCell][iN] = -cf/l1**2.\n",
    "    #M[iCell][iS] = -ctS\n",
    "    M[iCell][iW] = -cfW         \n",
    "    M[iCell][iE] = cfE\n",
    "\n",
    "    # Fracture tip 2\n",
    "    i = int(n/2.)\n",
    "    j = int(60.*n/101.)\n",
    "\n",
    "    iCell = i*n + j\n",
    "    for k in range(0,nCell):\n",
    "        M[iCell][k] = 0. \n",
    "\n",
    "    iW = (i-1)*n + j\n",
    "    iE = (i+1)*n + j\n",
    "    iS = iCell - 1\n",
    "    iN = iCell + 1\n",
    "\n",
    "    l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "    l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "\n",
    "    cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "    # Removing ctN allow a best fit for high matrix permeability, steady state regime\n",
    "\n",
    "    #ctN  = coeffff* (  -phi*l2/e*(cyy[i][j+1] + cyy[i][j]) /2. /(y[i][j+1] - y[i][j])  )\n",
    "\n",
    "    M[iCell][iCell] = s + cf/l1**2. + cfW - cfE #- ctN\n",
    "\n",
    "\n",
    "\n",
    "    #M[iCell][iN] = ctN\n",
    "    M[iCell][iS] = -cf/l1**2.\n",
    "    M[iCell][iW] = -cfW         \n",
    "    M[iCell][iE] = cfE\n",
    "\n",
    "    # Pressure boundary condition on top and bottom sides  \n",
    "\n",
    "    for i in range(0,n):\n",
    "        for j in range(0,nCell):\n",
    "            #M[i][j] = 0\n",
    "            M[(i+1)*n-1][j] = 0. \n",
    "            #M[nCell-1-i][j] = 0\n",
    "            M[i*n][j] = 0.\n",
    "\n",
    "        #M[i][i] = 1\n",
    "        M[(i+1)*n-1][(i+1)*n-1] = 1. \n",
    "        #M[nCell-1-i][nCell-1-i] = 1\n",
    "        M[i*n][i*n] = 1. \n",
    "\n",
    "    # Flow rate boundary condition on left and right sides    \n",
    "    for i in range(1,n-1): # 1 to n-1 to avoid the boundaries where pressure is imposed, this allows avoiding an issue related to 2x pressure at the first cell\n",
    "        for j in range(0,nCell):\n",
    "            M[i][j] = 0.\n",
    "            #M[(i+1)*n-1][j] = 0 \n",
    "            M[nCell-1-i][j] = 0.\n",
    "            #M[i*n][j] = 0\n",
    "\n",
    "        M[i][i] = 1.\n",
    "        M[i][i+n] = -1.\n",
    "\n",
    "        #M[(i+1)*n-1][(i+1)*n-1] = 1 \n",
    "        #M[(i+1)*n-1][(i+1)*n-2] = -1\n",
    "\n",
    "        M[nCell-1-i][nCell-1-i] = 1.\n",
    "        M[nCell-1-i][nCell-1-i-n] = -1.\n",
    "\n",
    "        #M[i*n][i*n] = 1 \n",
    "        #M[i*n][i*n+1] = -1  \n",
    "\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:30.415289Z",
     "iopub.status.busy": "2023-10-06T07:42:30.414588Z",
     "iopub.status.idle": "2023-10-06T07:42:43.118807Z",
     "shell.execute_reply": "2023-10-06T07:42:43.118132Z",
     "shell.execute_reply.started": "2023-10-06T07:34:33.533416Z"
    },
    "papermill": {
     "duration": 12.7153,
     "end_time": "2023-10-06T07:42:43.118970",
     "exception": false,
     "start_time": "2023-10-06T07:42:30.403670",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "p0 = 0.\n",
    "pTop = 5\n",
    "pBot = -5\n",
    "L = 10.\n",
    "e = 1e-5\n",
    "\n",
    "\n",
    "n = 101 # Highly mesh sensible for large time transient regime, this may be avoided by removing time dependent term in the fracture\n",
    "#x, y = np.mgrid[0:L:n*1j, 0:L:n*1j]\n",
    "x, y = np.mgrid[-L/2.:L/2.:n*1j, -L/2.:L/2.:n*1j]\n",
    "\n",
    "t = 1e6\n",
    "\n",
    "muf = 1e-3\n",
    "\n",
    "c_fluid = 5e-10\n",
    "kappa = 1e-18/muf\n",
    "phi = 0.1\n",
    "\n",
    "# Biot's modulus\n",
    "M = 1./(phi*c_fluid)\n",
    "\n",
    "# Diffusion coefficient\n",
    "c = kappa * M\n",
    "\n",
    "\n",
    "\n",
    "cf = e**3./12./muf/c_fluid  \n",
    "\n",
    "cxx = np.r_[c*np.ones((int(n/2.),n)), c*np.ones((n-int(n/2.),n)) ]\n",
    "cyy = cxx\n",
    "\n",
    "# Force\n",
    "nCell = n*n\n",
    "#f = np.zeros((nCell,1))\n",
    "f = np.zeros(nCell)\n",
    "\n",
    "# Pressure boundaries condition \n",
    "    \n",
    "for i in range(0,n):\n",
    "\n",
    "    #f[i] = 0\n",
    "    f[(i+1)*n-1] = pTop\n",
    "    #f[nCell-1-i] = 0\n",
    "    f[i*n] = pBot\n",
    "  \n",
    "f = f.reshape(len(f),1)\n",
    "\n",
    "# Solving the linear equation system by the fast Laplace-Stehfest method\n",
    "N = 2\n",
    "sum_p = np.zeros(nCell)\n",
    "\n",
    "for j in range(1,2*N+1):\n",
    "\ts = j * np.log(2.) / t\n",
    "\n",
    "\tM = rigidity(n,s,x,y,cxx,cyy,cf,phi,e)\n",
    "\t\n",
    "\tp = solve(M,f/s)\n",
    "\n",
    "\tsum_p += Vfunction(j,N) * p\n",
    "\n",
    "p = sum_p * np.log(2.) / t\n",
    "\n",
    "p+=p0\n",
    "\n",
    "p = p.reshape(n,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:43.140513Z",
     "iopub.status.busy": "2023-10-06T07:42:43.139765Z",
     "iopub.status.idle": "2023-10-06T07:42:43.510505Z",
     "shell.execute_reply": "2023-10-06T07:42:43.511005Z",
     "shell.execute_reply.started": "2023-10-06T07:42:02.095565Z"
    },
    "papermill": {
     "duration": 0.385139,
     "end_time": "2023-10-06T07:42:43.511174",
     "exception": false,
     "start_time": "2023-10-06T07:42:43.126035",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEOCAYAAABIESrBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dd5wlVZX4v6dfTw/DMEOaUYIkAYEBBXFWBSUMLAIuSxQRAYVdCSZUTAvKijmCoKvguKaVFYkuwo+MIigoDkEEAwgSHBCZYYAJPaG7z++Peq+7+nW9rnTrVd165/v59Kffq3DvrXpV99wT7rmiqhiGYRhGFvrKboBhGIbhLyZEDMMwjMyYEDEMwzAyY0LEMAzDyIwJEcMwDCMz/WU3oAjWFtH1ym6EYRhe8BQsUtXZWc/fRkRXJK/relU9IGtdVaSWQmQ94OSyG2EYhhecBY/lOX8Fyfubs2BWnrqqSC2FiGEYRrcQersjreW19wHTym6EYRg9Qa/3N7UUIoZhGN1CgCllN6JEailEGsCMshthGEZPYOYswzAMIzOmidQQH22UvfwQGvViTdkN6DKmiRiGYRiZMU3EcMpM4IWyG+EJ04DBshvhCXH3yu7lGHHvoOt31EfLh0tqKUQaBA9KWZRZt2+0Xr46jeSKMufEdVRV6ciq0I64d9DlO2qaiGEY3mMacHKKuFe93JHW8tr7MG3AKI+yRuL2zCfHNBF31FKIGEavMQNYWnYjehSLzqohdXF01fLHMQqjjAm2QyXUWTXq0t9kxfopwzCMHJg5q4ZY2pNi6eUXpix6bQKfTxRpzhKRM4DPAt9Q1fcUVE0uailEDMMwukVRmoiIvBY4EbivgOKd4Y0QSSORy54nYrinStpPGVpBL9vcq04RmoiIrAv8L/DvwH86Lt4pXqyx7otENgyj92hpIkn+UjAfuExVf+awqYVQeU0ki0QW6jFym1L5X6c3qcOz5Yo1Fp6Vtr+ZJSILQt/nq+r8ceWJnAhsAxznon1F40M3NSqRRaTSap1hGL2HkGLAN8QiVZ3bsSyR7YDPAXuo6moX7SuaSguRNBJZRE4CTgJ4ScHtmoyZ0+GF5SU2wOhJpk2FwVVlt8IPXL+jItCfXIjEsRswC7hfRFrbGsCeInIKMF1VK/VLV1aIpJXITZVwPsCuDdGZJdocZk4vr24XJH4hjMowNBQIkm5TRp0uGPeO5hQoIjClka+MEP8HLGjb9j3gIYL+sHLaSZW7C+8ksmGUxbTpMGgacCJc36tUmkgMqvoc8Nz48mU58Kyq3u+mFrdUWYhklsh9fcGDYpRMlZ+uquDQMW3PfHKmudZEPNXIXFDZ19xHiWwYpbEOsKzsRvQoBWdgVNW9iys9P5UVIrnoA9YquxGOqecvlRx3Nuf8DJfdgA6s06V6LKx3PD2exterS6+6RDYMo0fxqid1Sz0vvQ8w+3A89fz1jTBFaA09bP+PRKiWphyDiOwAHA3sBWxJMFfyGeBu4Frg8jRBS16kPTEMw6gsLXNWkr8SEZFdReQm4B7gdcDtwFeAM4AfAEqQn/BJEfmoiCQaLtRzLCoUN1qq5x3Lht2L4smrSdRJa6iqL6bI/sYtPwG+BBypqks6HSQiuwEfAD5IEAk7KdYNGIZh5MEfx/q2CSdu3wHcISIDSQr149LTUnefiI+/mo9tLoo0I+oqjXDL1gSqdC/CeCJE4gSIiExR1TVJj29hPhHDMIy8NBL+VQQROVVEjgh9/w4wKCJ/bqacSowJEcOoA3XWvF3jej6NJ471Nk4liMhCRPYE3gy8FbgXODtNQdW6LFc06N7EK5+p0MgoEh+eziwmnqLMMlGCpKoTI8vGZf/giTmrjU2BR5uf/xW4VFUvEZHfA7elKcg0EaO6lG2D94m6ZWgoEtdCvBWdleSvOrwAzG5+3g+4ufl5DSmfJv/kZxKqGHLn+52uutYSJsu9ziqwynjOOmkXvgiSKgwOXP5ufmoiNwDfFpF7CNZsura5fUfgr2kKMk3EMAwjD376RN4N/IpguY03qeqzze27AhelKahal+WKBr3haKzqr+eT1tKibM21qr4LS5sSj2dpTwBU9QXgvRHbP5G2LNNEDMMw8uChJiIic8KhvCKyn4hcKCKni0gqkVihy3JI3ScbuqKsX7+eT93klOUHsLQpxVNFH2w83wHOA/4sIi8BrgRuITBzzQROT1qQaSKGYRh58FATAXYgyNoLcCTwG1V9I3AcQYbfxFTrslwh+BOp0qKKNlXTVJITN+IvaqRaVr2TUVX/TlH4GZ3VYGyJ8X2Ba5qfHwZenKYg00QMwzDy4p8mcj/wThHZg0CIXNfcvimwKE1BvSdE4nwlRe/3ibhRrN2rMYq+F3l/iypR9jvo+l61orM8yp0FfBQ4kcAPcpGq/r65/WDgzjQFVUs2uiIuxNe3hzQtLh/WqgiSKj2pnUxIRdyLcF1Rv0V4W97nrptmqLLfQZfvqIfmLFW9VURmAzPb1hb5FrAiTVmeXbphGJGsAywruxGeMBVIvPhrAvyMzkJVh0WkISKvAe5V1VWq+mjacuopRMr4UX28k0W3uRvqe9HXEKV1uHy2XI7+o5IKViHFSFq60eYeT3siIjOA7wJHECyLuy3wiIhcAPxdVc9KWlbv+UQMo454OBKuDX6G+H4R2IQgzclgaPvVwGFpCqrWZbmijpMNy/iliqqzDAdj+FrKWFmwKH9DGp9JEXUWgW8C0UNNhMCBfpiq3isiGtr+R+ClaQry79INwzCqRrUir5KwPrA4YvsMUg556ilE6qiJZKHoX9dl+d18CYse6brUOoqeTFilRbV8xU9N5LcE2si5ze8tbeRk4PY0Bfl36YZhGFXCz+isM4DrRWRHAjlwWvPzq4E90xRUTyHShx9pT6qkArt8EvKW5eNTmddP0On8LJ1TGT6fTvRCChQPNRFVvV1Edgc+RJDqZF+CXFq7hSYeJsKzSzcMw6gYHgoRgKaweHvecjy89AT0ik/Ep3keWdrqw9PZGvWnGdG7nHtS9Ejfx3km3caTRalEZIOkx4ZWOozFh9fUMAyjuvijiSxizIHeCWkek1gs+nHpPjEdWF52Izwh7l4NYU9oi7WAlZPst+dujLi0JkWkPfHBBwvziii0nq9oH+VGS1TdlObqV8+qwofr9z0BY9LnLMs5MN5cFddRxd3LqofzujSdxbXbddoTR+YsETkdOBzYjkDU/Ro4XVXvz1u2qv4ibxlRWNoTw6gD/oWYlkcRqeDdpT3ZG/gmsDuwD4FovSmNPyNRk0WGReRFEds3FJF6TDbMJZHr6Fjv5i+Vpa4sI7Gs11TGU5tVk8jieM/iLB9KUIcrQVOUs72bgtB1/+DomVTV/cPfReQ44HngdcBVbmoJiu6wfSpjKx4morJChDGJ/FuCC/4UgUSekyZywDB6AvOJlEex0VkzCIbFS+IOTIKInNb8qMApIhJeQKAB7AH8KU2ZlRUiuSRynsmGHoTqJabsCYRpzklz38vWROJIo0lkGdVnXRTLVf1VNZ2VNbExXXTWLBFZEPo+X1XnT3L8ecC9wB3ZGjeB9zb/C/AOxt+11cCjwClpCqysEInAqUQ2eov3fy74f+4Z5bbDqCHp0p4sUtW5iYoVOQd4PfB6VXUiIlV1q2bZPwcOb1vVMBM+CZFJJbKInAScBLD5RvjtE+nWr5JX60rTTpdaSYb7c+9DzQ9Zn4u4EXtW7SOu88k7MbGbo/NenZhYwDwREfkq8BZgnqo+4rZ0UFVn4b5eCJEkErmpEs4HmDtH4ibUGIZhuMGxEBGR8wgEyN6qmso/EVPu1wiCk5Y3P3dEVU9NWm7lhUgmiVzH6Kw4ivglXcwDKWC/ZmlX8xzN+FzEBj3GjcLjtI+sCRiL0DSyahRV9ZUUjUMhIiLfAI4DDgWWiMhGzV3LVHVZ5zMT8XJgSuhzJ1INwistRIqSyIZhGE5xF5Dzrub/m9u2fxI4K0/BYRNWT5izcklkX/L7l3X3u5SqPU5jGI4pp9P+ocbkc2SH+yeeONxYA8DSdadM2AfQGJp8+N0/PNLhvOaHDs9b3P7E07pc+BuSluHy3ekFP4lDTURVO83fqCyVFSIUKJENwzCc4cugtQ0ROYpgHZEX0Za9RFUPTlpOZYWIjxLZcEtjKF5bMQJ0KojLpIJ1Zh0gr3chjD9ZfEcRkS8D7wd+DjxJSj9IGM8uPSENggel5mRyMCcgS8cdd042E1T0BQ7HGKCHIvYPN+1GKxrTIs/pb0xuV2p08GA3hiduH2caixihdjKNjZ0fsTFmpNsYCgTJpGQYLafLouQRLvsHD4UI8DbgaFW9LG9BloDRqCxRHbQRzYrp9ionZci16cltAsZu0Ucw7y431bosR6gU8KDkpEpmmTitIIooTaHjsR00iNH9EZpCpPbQ4fGMOj9OOxlmEIBlzOioVbSI2t9RE4lQG6K0mvD5UdntxgnMWK0j2lsdJ0jiNKBuEalpdRnX/UNRVoECmQ8ciwP/coW6NsMYT4OhjoLEGM/qqVMZWGVOkTLQPljtx6JUYdYD3ioi+wH3AWvCO2s12TALIw1hxfToUM46ETfiz1xuwqD3KO2hc5kRPo8O58dpGnFaS9T5IzwHwArWjqwzrCk0IuJS+8ftD2saAxO2hTWNaK0mYije4Vb2R2lAjYGI04dZPXXy4XWq/N6tcoswKVbMSpDtzoyhkka7r4Y2CMxhzJy1fdu++kw2NAzDqDoqksLcm09guaInJhvmQRFWR4zWukGa0Xk3yWsWSqqdJDk2zicSp7VEaSVx/pPWcUuZMW5/nKbQ0gQaHTWRzud0ridak4gi0n+SM9dJpCbU+eCuEal1dYX8HXtRVgEfqKUQMQzD6BaKpBpkVQER+elk+2sx2TAPI/R1tH3XgaIfWFc+kTjtp5N2EVV+J5/H8KgmEu0zaZ031Pz/HOtF+j86+TwakZpIlNYxMOGc9s9RxGkFSUfn6bSTbFp6Xg0ojvIMPflmHipSWQvEJCxu+z4F2BnYDLgiTUG1FCJGPehn2MeX0+gxFGF19aIFJkVVT4jaLiJnA0vTlFVLIRL8qMX4RHzt1KrkE0kaXRXnBwmfF96/KvTbt17uNc0M2IvYkKmhMe8AQVhseJtLTaVFZ41iYnRXmNJ9IhWiPJ/J5PhozpqEbwG/JMX8kVoKEcMwjG5SIyGyXdoTTIgYlaXBcJ1ezkJZzQADFQkfrTorWJu1WeGsPB99IhErGwqwMXAg8N00ZdVSiJhjPR9ZXog0TvSobVkd6y3TVXj/CsaSLA42n4OVBFOKH2UrXszTo/unNTuTVloUGG/uaZm5hmMd64T2TzR3xU1AHH9+enNTP8MpTLjVdKzH4VJEuuwfAnOWd11p+8qGI8AzwAcwIWLUBXOsJ8f16NpITpE+2KKwyYYxjCDjnKud8HD04IS8mkzS8zsdlyVtyXjtY+qE7VHOdIgecYa3RY2uw072Vrumsip0zlhbVo9uy+ZYjx7dJ5+MGCZudF22JlEG3QgiUPwNuHGB5Y82DMPIRWDOSvJXaitF/ltEtkhwnIjIMSJyXJJyazkUV/pGbeF1pFvO5qyjqzTJFqP2R6UyiQv3DWsfYdNCnFYTZ4ZoaRDhc+JDgCfWmSkpYwK6GfbqpyZTvJnJoxDfvwH3ichvgJ8CC4CngJXA+gRJGV8PHAU8CpycpNBaChGjHlgqeMMXfBAiqnqWiHwDOJFAQJzXdshS4CbgBFW9IWm5tXxD6x6d5Ypu+UYg26JTcRMMw8eENYqhDlpL1LaxY6NHrK2orE4RWY3RtnbyeST1j0TXHzf6j4tY8lN78AuPNBFU9Rngc8DnRGR9YHNgGrAIeFhVU6+1nkqIiEgfsJaqWhiIUTg2T8TwAUXGBXv4gqouAZbkLSdWiIjIgcDRwF7ApsEmWQncDVwHfE9Vn8zbEJeM0OddyF1euhkdktTElCUlPCRP7w5jGkhcKvnObYg6diD0aeJYP3xOlE8kTBb/SHQ53UtVUtX0IlXFJ02kCDq+bSJyKPAlYAZwDYEK9CQwCGwA7AT8M3CmiHwfOLOpKhmGYfQMJkQ6czpwGnCNqkat6XgJgIhsCrwPeBtwtvMWdpm4SVt599eJuFQbqxgYN+di4v6p4+ZftNNLjvW4eznI2qOz66PPnzqaTLLuVPEd7eV5Ih3fUFV9TZICVHUh8BFnLXJAXsd63Ll591edNKOqKLNh+IWKm/QZZUtOkyIlLhw46uUOb4syN41nYpbd8OexIIDkjvUokqQtmSxsvcGws3TkPmTgqtI76mnaE2fYZEPDqAG+D1y6SZJsFmlombOS/FUNEZklIq8RkcwjkMTiU0QOA+YBL6JN+Kjqm7M2oAg0YdqTutPN0VGWF6TTOXHJGMcfO3E9EZdrn0TVH7+Genh7svKzONvbySJIfF1bJC8u+wcf+xsRmQF8B3gTQeaWbYFHROQC4O+qelbSshJpIs3Vri5mLPPjcNufYRgl0mvRiFVCPUl70sYXCaJtd4VQCmu4GjgsTUFJr+rtwJGqemWawstipOZpT6pCXvU86xrtyScQumhDsH+8H2RiiG+ncvKGy6bRSrIJEhM+LqiiqSqGg4HDVPVeEQlPMPwj8NI0BSUVIiuAP6Up2DDyUsXorEfv+Du3fvV3/PWXT7F80UqmzpjCS3adxdxjt2XusdvQ1+iOm/Gd8n0O+NgrOOQzu3alPqMznob4rg8sjtg+g5TWpaRv6BeAj4jIyapaeSOqjzbKKtGtjjvri5fFZ9JpW5wmEea2c+/l6tNuY+t9XsKBX3wd628xg8Elq3jwhse59J2/ZMp609nxkPGDuE5RWy4YYsqoH6RXfRtVwFMh8lsCbeTc5veWNnIycHuagpL2Ft8GDgIWisiDwJrwTlXdJ02lhuEbj9y6kKtPu43d3/MKDvnaXuP27XjIS9nztFeyevmayHOHVg3TP9W7TsZISBFpT0TkXcCHCZasfQB4v6re5rCKM4DrRWRHAjlwWvPzq4E90xSUVIhcAOxBkObkacakVuFkuZnmE/GbvBO3sqQ9CY/ko+q/+Qv3MG2Dtdj3S/NYFVH+eltvAMCd33+Ay064kZN+8SZu//rv+MuNj7H+ljM59d5jeeK3f+fWLy7g8V//nRWLB1lv85nsdMQ27HPma5gybXyZD/zkL9z2lbt46neLkD5h9vbrM+/jr2bOwVuH2t4Y9YOsXiFcdNQ1LFzwNCdcdxgb7zw78T0w8uFaExGRowgy7L4L+GXz/7UiMkdVH3dRh6reLiK7EfStDwP7EqSy2k1Vf5+mrKRv25uBw1X1xlQtzUk3bqZhxDEyPMKjtzzB9oduw5S1kr0yFx9zHTsfvR3HXHYQw0NBwofnH1/KxrvMZtfj5zB1xgBPP7CYn33qNzz7yPMc/eM3jp57+9fv5apTb2HOoVtz5A/ewMA6U3jy7n+w5NEXIuta8exKfnDQlSxfNMgptx/FBlutm/+ijVQ4NmedBnxfVb/d/P5eETkAeCdBJpFciMgU4ELgDFV9e97ykgqRRcDCvJVlINPNHKHPy6yaRnlMpr0sW7ScNYNDrLvFepNEjAXb1zAFgB3etD37fGnf0f2rgG2P2Iltm99VlY1etxWNmdP5yduu5oBvHMDaG05j5QuruP6M29n+sJfx5ivGIi232H+70XJaDNHg6cdXceH+lzCwzgD/9qvjmD577R5JflIdFHGW9kREBoBXAV9p23UDsLuLOlR1jYi8AQcCCZILkU8AnxKR41V1mYuK4+jGzTSMotj+sJdN2LbyhVXc9tk7+MNlf+b5J15gZM1YSrrFDz3L2htuyhO3L2T1stW86qSdY+t45g+L+c7uP2T2nFkcdcVhTF3HgknKIGXak1kisiD0fb6qzg/vBxoEboMwTxMkvHXFFcDhTOxfU5P0yj8MbAk8LSKPM9Gx/oq8DYkg1c0UkZOAkwDW3XxmAc1JxiDTmDZu7o7Ribikg1Vh2obT6J/Wz3OPPZ/4nBkbT5+w7coTruGRmx5l3qf2YKNdXsSU6VNYeOdTXPPuGxlaGURxDS4Onp2ZL5kx7tyoe/XYrU8wuHiQ/c/eh5F11gV77oD4d7CIdzSFOWuRqs5NcFy731kituXhceDjIrIHwTK5y8dVrnpO0oKSCpHLkrfNOYluZlOazwfYeO4mOsi0LjQtmjLrdkE3V8PLMkEuTfvijk308vc32GLvLXn4xsdYsaqP/qmdX5tRs5YMjF9tceUQf7ryIfY4ay92fd+YIv3k74M1gYboDwTFrGAAtGThStbfafy9ab9Xrzx5V1Y9v4orjr2ag/vXYrsj5sRfS48Q9w66fEcVcZkxYBHBPI2N2ra/iIkD6jwcT7Ag1Suaf2EUcCtEVPWTSQt0SLdupmHEstt/7M6Fe/8PP/vwTbzhawdM2P/cX5ewemlnrWpo1RA6rDSmjJ+MeN/3fzfu+6a7v4SBdQa4Z/7dvHT/rZkUEfb/rwPp6+/jp2+5jH/90RFsf+SOyS+qR4lbgiAtLn0iqrpaRO4C9gMuDe3aD7jcSSVBPVu5KmuyRakkzXq7aY+PI8/NNMd693E52S0uVUhYe4jSNKJStXc6f7Jt4e2b7vlS9jlnf3522vU888fFvPz4XZi5+bqsXDLIYzf/ld/99938y4+OZM3omu9Txk14lXUH2Pi1m/Hrs3/D1I3XY9qs6dz/3btYujBwMa5pHT9jgN0/vz+3vPcqLjniCrY/ZhemzJjKM/c+Rf9a/ezy3jEtpjXZ8HXnHsZwY4Cr3no5q0f62e6oIqzL9cJl/1BAKvhzgB+KyJ3Ar4BTgE0IplpUjsmu/E8i8mngMlVd2ekgEdmBIIrqEeDzjtvn1c006s0/vX83Nn71piz46q/5+YduYHDRCgZmDLDR3E3Y71sHs/W/bscD/3Nvx/P/5aIjuemdV3Hzu6+mf9oUtnvzTux03q785KALxx23y3t2Z/pGM1jw5Vu59piLaUxpsP4Os3nNmZ3n9L727MOQ/j6uO+ZidGSE7Y/exdl1G/G4DPFV1YtFZEPg4wTz4+4H3qiqj7mqQ0S+FtOGUxOX1Ul5EJG9CJbH3Q64icD58hSwkiDvyhzg9cDLgK8Bn1fVpUkrTtzAYLLhRxi7mR9Q1VsnO2f23M318AUfdN2USlPFdbFd+i46aTpxadlbn/9n7x8CcMItb41MGx+33nuYuFFnEanou0mvrdI3X95/V0JndyQvmruZHrngA4mO/aZ8MFddrhCRn7dtmgJsT6BY3J0mC8lkKxv+AniNiOwOHA28BdgCmEbgr7gH+C5woao+l+oKUqCq3wS+WVT5hmEYeXDpE+kWqjqvfZuIrEWwxkiq9CqxhjxVvZ2UCbnKJoiW6C2fiMtAWVfRWWn8JHEaRZwfxKUmNjROU5m46FWYNMkgo893Y0uvokbTK9Slv1HVlSLyWeB6UrgMqpVn2zBCDNPoarixYWTB0yy+nZgNrJPmhFoKkRH6WOH5XI12quTzSN6xj0Unxc/XiNJEhkZH9S6vP8onEqV9hPd3WvQqOtV8soWu0lClTso30003qNLvkwQROa19E4Hf+RjgmjRl1VKIGPVgmH5bJ8OoPD76RID3tn0fAZ4BvkfKKNtaCpG62CjDFJEcxKWpKK6zj9MkoiKqIKwJhPeP1dUqt9N84bhrHBotf+xVGDfTPGJ/3DyTuFFp1mWBs+DbCNlHCpgnUjhdmWxoGIZhxOM47UlpiMg2wN8mmxcYRaIFoUXk/0TkIBHpzgLShmGkom4+wCJZ4XjBupY5K8lfVRCRz4nI25ufRURuAh4EnhKR16QpK6kmshy4GHheRL4PfE9VH0pTUTcJHOv+r2xYtj8gnTN7cid62BwXnapkoolq/P6wOSv8MrZyII3VHx4TjpnBpPm90WZ66m+2b6IJK9g+dfS8sXMmD+vtlrO9vd4yBIlvZpwWrvsHD+/DMcBRzc8HAjsDr21u/wIwYR5JJxJpFqp6DIHn/tMEadj/LCK3isjbRMSGQIZRMnXzARaJa9NTK8Q3yV+FeDHwt+bnNwKXqOqdwNeBV6YpKLH4VNUXgPOB85sLur8D+BbwdRH5MXCuqv4xTeVFMUJfLWyUnd3F6cniRO/kzE9aVrpUJeHR/8R64jSV8Jp/4bKiTAhDEVpFlPYRfB6IOCe5k32ydrSXNRlJOqCkgqRinVkpuOwfPJ0nspggA8nfgPAqh/201PaEpPZxiMgmwCHAQcAQwVojmwH3iciH0pZnGJ3w8MUsjTqYb33GN58IQSb0H4nIjcAGwHXN7bsAf0lTUKJhUHNh90OAfyNIxX4PQXLGi1rL5YrImwkWhcq93GJe6uMTKWeCYRZfTLR2EO0nidI0ourvZ3hUkDTGaQL5EjuGR/+tdO1R2kewf3KfSFQI8Pi6kqedby8zLWU88yboW5YP78yJpwGPAZsDH1HV1sqGGxNYnBKT1Jz1FIGK8yPgP1T1vohjbiRYKcswnDBEo1Iz9Q2jE74JU1UdAs6O2P7VtGUlFSIfAC6dLH5YVZcAziaw5KEucdt5ya7JTLx3cWXFRV+NL2tM04je39I+hkdH+J18IlE+jTBRdYSfjdYIMhzZFB5VZkmLEiZOU2kvJwkuOyzfOr8q4qNPpLnUx0pV/U3z+/EEfu4HgA+2LExJSBqd9cO0E1AMIy++vZhGb6J46RM5l+bS4yKyHUGQ1H3AbsCX0xTkXXBzEswnkqUuV36QcJmTaxpZ546Ejx1oRmVFzf0I6gi0jpVN7eIJNov0eaSZJ5JlTkiaVPJhsnQ83ZyzYIIe8DDtCbA18Pvm5yOAG1X1Xc2JhpcD70xakHdXbhiGUSV8NGcRKFCtRu8L/KT5+e/AhmkKqqUQUWQ06qaXyTrjPIq4Wehxx0ZpOp19IsPN/xPnjrSX3/qdp4aOWBXyabR8HUNMAWAxG8YuOhV+dpIuSpV3ed2yfCIVM7F4iaf9zW+BM5shvnsAJzW3b0kgSBJTSyFiGIbRLXzM4gu8nyDa9hDgs6r6cHP7kaRcyda7Kzd6B1vZMDmrmMrU0Kx9ozODrM00Vjgt0zdzlqreD7wiYteHIN1LV0shMkIfg+ZYT7UGSZa60q0hMsP2pgcAABNlSURBVHHiYd6kjOPLmjphW7iMNc1HfQnrjTs/ylxVRLLFNOauqPYlIWpiZNH41nm2cNk/eOoTAUBE5hI42a9uTjhsQLoom1oKEaMe2MqGyRkcmca0vsGym9GTKMLwiF9CREReDPwU+CcCJ/u2wCPAOcBK4H1Jy6qlEBkZ6WPFivwjjUa/3x1Yf7/DlQv7Jk9rMu7YhJrG+HM6Od4HJj1/zAk/ebLH1khxGTPG7Y9zfMdpFUkd553On2xbx/M7dFiDI+kTag8N+dX5tTM8VH4XpiPCqpXepT35KmORWI+Htl9KkMk3MeX/AobRAUt7YviAqjDsnzDeF9hXVZeIjEva+zBBPq3E1FKI6IiweqWLkLvqhO01MmgVqXwiGcqP09TiNKHx2k2+ZI3j6o1Muhi85EvbNJEWaVK1p5pMGGPmiNMEsoy0i+7QPOwwi0W9vCfTiO4iZhOYsxJjy90ahmHkQFUYWtNI9FchbgWOD31XEWkAHwVuTlNQTTWRPlYv63J0lkP/QxH0NfL5d/qnpPCJ9E/0iUTt70SUhtNJq4n21YT2N7WSkaam0J4OJ3bZ2ghNIkp7iNMY4kaqaUayeTujkeGKv/bejeql+vd0Ih8BfiEi/wRMJcjouyOwLvC6NAV5d+WGYRiVQvFO8KnqH0Tk5QQ5slYBaxE41b+hqk+lKaueQmREYGW3L62L9fVr6lNG0vh3Ikb9q+PmscVoF3Ga0HhNJyIVfSdNJKEGONTUKJauiPaJdCJKQ4jTGuI0hdhRa5YOqRud2FCqVVN7h1L6m+w0Fxn8JfA2Vf1E3vL8uXLDMIyq4tFsAFVdIyJbEehQuTHHelrilmpJvJRLDxATITeyyrvY+uJYFnMvlsWM91b2kJaQ9x10/Y4GC4ok+6sOPwBOdFFQPTWREYrtzLM8pE7vdM4OI7YtMQdE7o7Y2K+jgiTSnNY/PCpIIs1lKYIV4sxlI0PBeGnF0vgJeamcpHFmpNj9od+yaEFSlU7MRTuqJEhaQsQvpgPHiMh+wF3A8vBOVT01aUH1FCJGPRhqVD7qrTIshQ5TYIx2VhK4kV2hwBqH5SVERDYAPgnsB2wBLAKuBj6uqotjTt8BuLv5+aVt+1KZuSopRHLenEAT6bXFfIvwq07JeN7oUxUzMu5PovEke0RH+mMCB5qO9ZHlCUK/0ziQ40agSUeoWUey4c5racYy2ukFue2yf1AoKYHyJsCmBOG6f2h+/iZwEfCGyU5U1XmuGlFJIUKOm2PUiCGq+4RWjVWAuZjKoSRzVjOd++GhTX8RkQ8DV4vITFV9Ia4MEZlGkMUX4GFVTZ3Fs5KvaO6bU7RPpAjK+CXi6owN601RVxpNKYsG1B+jPbRG18sSaBlZOoQ05ow0I/00bYn6vcqw1fvnH8hHtXwiMwmehEkXTBGRqcAXgZMJYuoFWCUi84GPqmpiXa2SQqQDiW6OYRhGV0knRGaJyILQ9/mqOt9FM0RkPeDTwLdVNa5F5xNYdd4B3NHcthvweQLv2r8lrdcLIZLk5ojISbTWCZ69ee/5RMIU8atmKTNrO4qoa6tdgv9JNNQso8qsI9Fu1lVWuXUnnRBZpKpzJztARD4DfCymnHmqekvonOnAVcBCAjdAHEcCh6vqjaFtj4jIP4DLqaoQKfLmNKX5fADZZq6TSTRGjXj/uWW3wKgzbgXwucCFMceMrgEiIusA1zS/HpTQFLWcoE9tZyGQyi/SbU2kGzfHT59IFN1Mx5M1EiuKqmstaeiWpuCyE+pmuGkvRHLF4TgaVFUXEUSkxiIiM4BrCXwaB6hq0p7v68AnROT4ljO96WQ/kyovStWlm2MYhtE9SnKsN/vIGwj8xYcC05uWG4BnVXWyJYVeC+wFLBSR+5rbXk4gE6aLyE9bB6rqwZO1o5I+kZw3JxgddUPkVPLudSBvJFaWMuPIqmm51JqSknV0n3ekXkTn5JPvw4e2lhed9SoCYQDwYNu+ecAtk5y7iMD3EeavWRpR1W4wz80xDMPoHuXNE7mFjDmQVPUEV+2opBDJc3NKZxmwTtmN8ITlBBl8jHjiUnXYc5cc12lPwA+NqSAqKURyo5Qb4tst703Rv57L8juVtbzD9snOLyfFxBh5O4ys50c90+Gyinju6to5uk57Utf7lIB6ChHD6DWKGF3XFdda2wgpg2LrRT2FSLcc61mp4kqaZYf4FllOkbgagboopzW6LiGjbCKqFA7sOhV8la6ty/jwmhq9iiVgTI75RMrFzFk1o2yfSFVI8+uWEeLbiSpqanEUNRJN0zllGV33cOfnDE99IiJyIPBugvVE9lfVJ0TkHcBfVfXmpOXY8riGYRh58HB5XBE5BrgEeAjYijGDdoNkubdGqacmUnWfSBbKGJ0XPamvnk/feIruOMrwf/Sw/T8SPxfB+whwoqr+uKl9tPg18Kk0BfXCa2wYhlEsFdIyErItYyngwywjyBSSmHoKkRHcLRU6GWWk33BJEalQstbVS3Szw/GvcxtPVSPNwvjpE3kSeBnwWNv2PYGH0xRkr7ZhGEYeFD+E3XjmA18LmbI2E5E9gC8BZ6UpqJ5CZITuzGoue+Z0J8r4VcuKqCpaG+wln4N/o+lq4OE8EVX9koisC9xIME315wQ92ldU9RtpyqqnEDEMw+gWnpmzRKSfYGncc4DPAnMIInX/kGXJDRMihlEHlhKsjG3E43pipuJV2hNVHRKRK4DtVXUxsCDunMmopxDxcWXDqk+wKyOIIMnTWbZJsYwRaCcTWzeCSZLgg2nH0p78DtgGeDRvQfUUIkY9sLQnyVkFTC27EZ7gOlmlZ+asJmcBZ4vIJ4C7aMunrarPJi2onq+on5N/uk/RaVFc1u8r3excon4j/zq37mCp4P9f8/8VBFfQQprfE9tGeuE1Noz6YwkYy8PPEN95rgqqpxDx0SeShar/elX381QBl7Z0W5SqPDzziajqL1yVVfVuyDAMo9p4aj4XkRcTZPGdQ6BPPQCcr6pPpymnnkKkKgkY6zoS9z3di4/4Zy5Jhmcj+Eg8NGeJyOuA64CnGcuhdSxwmojsr6pRebUiqacQMQzD6BZ+hvh+BbgIOEVVRwBEpA+4ADgb2D1pQfUUIrYoVTKy/vplz80IU/QT7LtPwPf2+4J/93kX4PiWAAFQ1REROQe4J01BtiiVYRhGHjxclAp4nmAxqna2Ap5LU1C9hUiUX2RZyfs7HZtlf5y2Fbc/Tf1F3IvBDp+jiNsfN1s77/407Ys6tkrPXTefm7j9Zb+jLmg51pP8VYcfA98RkWNEZCsR2VJEjgW+TWDmSoyoavxRniEDc5WNcqWD8QPfjZF1DTwI45+tfDzVGj0XwxNyl6rOzXp6qv4mZ12uEJEB4MvAKYz1JGuA84GPqurqpGX53g0ZhgGWgDENRUzM9EzYNoXE+0TkdGBrgpnqf1HVFWnLqqcQKWKyYS+MmvNgYb/uyBouWpUEjEXjQrtznYDRkxBfEVmbQAM5lOCtvQk4VVUXZS2z3j4Rw+gVqhQx12u0QnyT/JXPJ4HjCXJn/RjYj8CElZn6aiLVcmJVn7xPgnVi5dP6DTwzrXiPXwkYDwf+XVV/DCAiFwK/EpGGqmYSc/UUIoZhGN1iBJ8WpdoMuK31RVXvFJEhYBPgiSwF1lOIeLbSWGqk7AYYPUn9AjndUQ1TVRIaQHvkVa6Ve+opRAzDMLqJPwJWgAtFJGyAXgv4toiMRmap6sFJC6y8EBERAa4F9geOVNXLSm5S+fjzwBqG0QVS9JM/iNh2YZ66Ky9EgA/ik7JoGIbRfRL1k6p6guuKKy1ERGQu8D7gVQQpixPil6erulT68TAAn8KCjGLI3k+6obK9hIjMIMjhcrKq/iPQ1gzDMKpGeYPWKvSTVZ5seAFwnapek+RgETlJRBaIyALIPPkyAS/k3G+M0StTrF2QN0NkL9Htd7Q1ZT3JH7Na/VTz76SclafqJ4ugq5qIiHwG+FjMYfMIYpl3BhInKlPV+cD8oJ6dtdjOvAxBUte8IiZIklMVQeJDjo9uvqOpZhsuikvAWGQ/WQTdNmedS3wkwOME0/LnAMva1LOLReQOVX19Mc0zDF95AZhZdiM8YRCY5rA858mzvOonuypEmkm+Ym1NIvIxguUbw/we+BBwZYKaqJ96X7frMdxjptTkuHyf3AqR7vWTbqikY11VFwILw9uakvYJVX2klEYZRqUxTaQ8yhm0VqWfrKQQyc8wNirrNnX12ZRB1lGtPfPl4FcGRtd4I0RU1WJ8DcOoINVZUKSMftIbIZKOYeoX9VP1n8p8Nr1D7466ozFNxDAMw8hMdTSRMqipELG0J4ZhdAvTRAzDMIzM9PagtaZCZASLVDEMozuYOcswDMPIhZmzDGfYpC+jDFyn8qgzrt9R00RqSNmTDS0Bo1EGloAxOa4TMNbhnmSjpkLEMHoN04CTU0QCRjNn1Yw6RkvU7XoM91gwSXJcvk917G+SU1MhYhi9xlJgRtmN6FHMnFVD6pj2xDDisGe+HMycZRiGYWTGNJEa0ts2SsMwuolpIoZhGEZmenvQKqpadhucIyLPAI91qbpZJFjK0jPqeE1Qz+uq4zVBd69rC1WdnfVkEbmOoL1JWKSqB2Stq4rUUoh0ExFZoKpzy26HS+p4TVDP66rjNUF9r6uO9JXdAMMwDMNfTIgYhmEYmTEhkp/5ZTegAOp4TVDP66rjNUF9r6t2mE/EMAzDyIxpIoZhGEZmTIgYhmEYmTEhYhiGYWTGhIhjJOA6EVEReVPZ7cmKiGwgIl8XkT+JyKCIPCEi54vIhmW3LS0i8i4R+auIrBSRu0Rkj7LblAcROV1EfisiL4jIMyJylYjsVHa7XCIiZzTfof8quy3G5JgQcc8HCdII+84mwKbAR4CXA8cCewIXldmotIjIUcB5wOeAVwK3A9eKyOalNiwfewPfBHYH9iFI3HSTiGxQZqNcISKvBU4E7iu7LUY8Fp3lEBGZC/wEeBXwNHCkql5WbqvcISJvBK4G1lNVL1ZAEpHfAPep6omhbQ8Bl6nq6eW1zB0isg7wPHCoql5VdnvyICLrAncTCJH/BO5X1feU2ypjMkwTcYSIzCAYpZ+sqv8ouz0FMRNYBawouyFJEJEBAoF+Q9uuGwhG8XVhBsG7vKTshjhgPoGA/1nZDTGSYULEHRcA16nqNWU3pAhEZD3g08C3VdWXvNezgAaBVhjmaWCj7jenMM4D7gXuKLsheRCRE4FtgDPLbouRHBMikyAin2k69yb721tEjgN2Bj5cdpvjSHpNbedMB64CFhL4SHyj3WYrEdu8RETOAV4PHKGq3vriRGQ7Ar/VMaq6uuz2GMkxn8gkiMgs4lM8P07g5HwbwcICLRrN73eo6uuLaWF6kl6Tqq5oHr8OcA1Bx3ugqi4ruInOaJqzVgBHq+qloe3fAHZS1b1Ka5wDROSrwFuAear6p7LbkwcROR74HuODUhoEwn4EmK6qq0pomhGDCREHiMimwPptm38PnAZcqaqPdL9V+Wn6ea4lECAHqKp3i3g3Heu/U9WTQtseBC732bEuIucRCJC9VfWPZbcnL01z6UvaNn8PeIhAQ3lArbOqJLayoQNUdSGBqWcUEQF4wnMBcgOBM/1QYHrTrAXwrEcmh3OAH4rIncCvgFMIwpcvKLVVOWhqUscR/C5LRKTl31nmk6YYRlWfA54LbxOR5QTP2v3ltMpIggkRoxOvAl7b/Pxg2755wC1dbU1GVPXi5gTJjwMbA/cDb1TVbq18WQTvav6/uW37J4GzutsUo9cxc5ZhGIaRGYvOMgzDMDJjQsQwDMPIjAkRwzAMIzMmRAzDMIzMmBAxDMMwMmNCxDAMw8iMCRGjtojIXiLyoIg0cpZzp4gc7qpdhlEnTIgYdebLwGcdJCb8NPBFEbH3xTDasJfCqCUisjuwPXBp3LEJuIZgzY4DHZRlGLXChIhRaURktog8JSL/Gdr2iuZ66ZOtYf9W4KZWNuLmeWeJyP0i8nYReVRElonI90RkoLkO+xMislhEzglrHU1N5hrg6CKu0TB8xnJnGZVGVZ9ppgm/SkRuJFh86SLgopilh/cALo7YviVwCHAQwRrylxMsUPV34A0E2sslBMkaLw+ddyfwsTzXYhh1xISIUXlU9XoR+Sbwv8AvgKnAe2NO2wJ4KmJ7AzhBVZ8H7heR64C9gE2bmYn/KCK/IkgyGRYiTwKbiki/Rys7GkbhmDnL8IWPAqsJFv86JkHK82nAyojtjzcFSIungQfbUts/Dbyo7bxBgnVV1krVasOoOSZEDF/YEtiMYKW7lyY4fhETFwoDWNP2XTtsa383NgBW+rpeh2EUhQkRo/KIyBQCU9ZPgQ8B54vI5jGn3QPMcdiMnYC7HZZnGLXAfCKGD3yawLz0z8DzwAEEqxXOU9WRDudcD/y7wzbsAVznsDzDqAWmiRiVRkT2Aj4IvE1Vn2uus308sAOBn6QTFwIvE5EdHbRhU2B3gjW/DcMIYSsbGrVFRL4AzFbVXBqJiHwZWFdVT3LTMsOoD6aJGHXmc8AjeXNnAf8AznTQHsOoHaaJGIZhGJkxTcQwDMPIjAkRwzAMIzMmRAzDMIzMmBAxDMMwMmNCxDAMw8iMCRHDMAwjM/8fVBRchlsOTecAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolor(x, y, p, cmap='jet')\n",
    "colorbar = plt.colorbar()\n",
    "colorbar.set_label('Pore pressure (unitless)')\n",
    "plt.xlabel('x (m)')\n",
    "plt.ylabel('y (m)')\n",
    "plt.plot([0,0],[-2,2],'k')\n",
    "plt.text(0.5, 0, 'Crack', ha='left', va='center', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:43.534132Z",
     "iopub.status.busy": "2023-10-06T07:42:43.533442Z",
     "iopub.status.idle": "2023-10-06T07:42:43.536467Z",
     "shell.execute_reply": "2023-10-06T07:42:43.536922Z",
     "shell.execute_reply.started": "2023-10-06T07:34:46.472643Z"
    },
    "papermill": {
     "duration": 0.017775,
     "end_time": "2023-10-06T07:42:43.537098",
     "exception": false,
     "start_time": "2023-10-06T07:42:43.519323",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def analytic(A, L, r, theta):\n",
    "    return A*np.sqrt(L)*np.sqrt((r**2-L**2)*(1.0+np.cos(theta)))\n",
    "    #return np.sqrt((r**2-L**2))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-10-06T07:42:43.572256Z",
     "iopub.status.busy": "2023-10-06T07:42:43.569323Z",
     "iopub.status.idle": "2023-10-06T07:42:43.738844Z",
     "shell.execute_reply": "2023-10-06T07:42:43.738314Z",
     "shell.execute_reply.started": "2023-10-06T07:40:05.309850Z"
    },
    "papermill": {
     "duration": 0.194098,
     "end_time": "2023-10-06T07:42:43.739048",
     "exception": false,
     "start_time": "2023-10-06T07:42:43.544950",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc8AAAFKCAYAAAB7BdRAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3hUZdrH8e8dCCUCNlAsJBFFARtCdEXXvrjYsKEgyNIRLFgX1Pha1kWsuKAooiKwRAUUXRfFtYsKliAWVtaCJohSRBCEhJY87x9nkpmElJNkJtN+n+uaiznPaXfOHOae0+7HnHOIiIiIfynRDkBERCTeKHmKiIjUkJKniIhIDSl5ioiI1JCSp4iISA0peYqIiNSQkqeIiEgN1VvyNLMrzOwLM9sYeC00s7OqmedwM3vXzArN7Cczu9XMrL5iFhERqUjDelzXCmA08C1e0u4PvGhmXZxzX5Sf2MxaAK8D84GjgUOAqcBm4IF6illERGQnFs0KQ2a2DrjJOfdYBeNGAPcAezvnCgNttwAjgP2dSiOJiEiUROWap5k1MLPeQDNgQSWTdQXeK0mcAf8B9gUyIxuhiIhI5erztC1mdjiwEGgCbALOd859WcnkrfFO9YZaHTLuhwqWPwwYBrDLLrt0ad++fTjCFhGRBLFo0aK1zrlWdV1OvSZP4GugE7AbcCEwzcxOds4tqWT68qdmrZJ2r9G5ycBkgKysLJebm1v3iEVEJGGYWX44llOvydM5tw34LjCYa2ZHA9cCgyuYfBXeEWaovQL/rkZERCRKov2cZwrQuJJxC4ETzKxJSFs34GcgL8JxiYiIVKo+n/O828xOMLPMwPObY4GTgZzA+LFm9mbILE8DBcBUMzvMzC4AbgTG6U5bERGJpvo8bdsamBH4dwPwBXCGc+4/gfH7AAeWTOyc22Bm3YCJQC6wHu/5znH1GLOIiMhO6i15OucG1HR84E7cEyMUkoiISK3U9922MWPjxo2sWbOG7du3RzsUqSepqanstddetGjRItqhiEicS8rkuXHjRlavXs1+++1H06ZNUbncxOeco7CwkJ9++glACVRE6iTad9tGxZo1a9hvv/1IS0tT4kwSZkZaWhr77bcfa9asiXY4IlLPcnJyyMzMBOgSjuUl5ZHn9u3badq0abTDkCho2rSpTtWLJJmcnByGDRtGQUFB2JaZlEeegI44k5Q+d5Hkk52dHdbECUmcPEVEJDksX7487MtU8kwgeXl5mBlV1fT1M02syM3NxczIy8vzPc+AAQM4++yzIxeUiMSd9PR0ANLCuEwlzzhhZlW+BgwY4Gs5bdq0YeXKlXTq1CmyAYuIxIgxY8ZwTJMmLArjMpPyhqF4tHLlytL3c+fOZejQoWXamjZtyvr166tdToMGDWjduny9fRGRBOUcfTdupHdREQ3CuFgdedZBya3PKSkpZGZmkpOTE7F1tW7duvS122677dS26667lk6bn59Pt27dSEtLo2PHjrz++uul48qftt2+fTsjR45k3333pXHjxrRp04Ybb7yx0jimTp1Ks2bNmDdvHu3btyctLY0ePXqwYcMGnnvuOdq1a8euu+5Kv379KCwM9mO+detWrrnmGvbee2+aNGnCsccey/vvv19m2a+++irt27enSZMmnHDCCXzzzTc7rX/BggWcdNJJpY+djBgxgo0bN9Zuo4pIYvvtN7joIrj8choE77IPywGokmctldz6nJ+fj3OO/Px8hg0bFtEE6ld2djYjR47k888/5+ijj6Z3795s2rSpwmknTJjACy+8wLPPPsu3337LzJkzOeSQQ6pc/tatW3nggQfIycnhzTffJDc3l549ezJt2jSef/55XnzxRebOncsjjzxSOs+oUaOYOXMmU6ZMYfHixRx++OF079699Oj5xx9/5LzzzqNbt2589tlnXHXVVYwaNarMer/88ktOP/10evToweeff86cOXP47LPPGDRoUB23mIgknA8/hE6d4Pnng21HHBG+5TvnEvLVpUsXV5mvvvqq0nF+ZWRkOLxOucu8MjIy6rzs6syePdt5H11ZP/zwgwPcpEmTSttWrFjhAPfee++VmeaTTz5xzjl31VVXuVNPPdUVFxf7WvdTTz3lAPe///2vtO366693KSkp7pdffilt69+/vzvrrLOcc85t2rTJpaamumnTppWO37Fjh2vbtq3Lzs52zjl30003uXbt2pWJ484773SA++GHH5xzzvXr188NGjSoTDyLFy92gFu9evVO661MOD5/EYlRRUXO3XOPcw0bOgfB1+WXO1dY6IBcF4YcoyPPWqrs1udI3BJdU0eE/Lrad999ASqtqjNgwAA+++wzDj74YK644gpefvlliouLq1x+48aNyxyd7r333rRu3ZqWLVuWaStZ57Jly9i+fTvHH3986fgGDRrQtWtXvvrqKwCWLl3KscceW+Y5zK5du5ZZ76JFi5gxYwbNmjUrfZUsc9myZVXGLCKJo9JLZmvWwJlnwujRsGOH17brrvDcczBxIjRpUvlCa0g3DNVSeno6+fn5FbZHW2pqaun7kmRUWULs3LkzeXl5vPrqq7z11lv079+fI488ktdff52UlIp/WzVsWHa3MbMy6yxpK1mnC3S/WlGBgpK2kmmqUlxczJAhQ7j22mt3GrfffvtVO7+IxL/y1YJKLpntvWQJf5o6FVatCk587LHwzDPgleULKx151tKYMWNISyv71FBaWhpjxoyJUkS117x5cy666CIeffRRXn75Zd566y2+++67sC3/oIMOolGjRmVuECoqKmLhwoV07NgRgI4dO/LRRx+VSaIffvhhmeV07tyZ//73vxx00EE7vVRuUSQ5lK8W1AC4qaCAU+++u2ziHD0a5s+PSOIEJc9a69u3L5MnTyYjIwMzIyMjg8mTJ9O3b99oh1Yj48aN45lnnmHp0qV89913PP3007Ro0YL9998/bOvYZZddGDFiBDfeeCOvvPIKS5cuZcSIEaxevZrLL78cgOHDh5OXl8c111zD119/zXPPPcekSZPKLGf06NF8/PHHDB8+nMWLF/Pdd98xd+5cLrvssrDFKiKxLfTS2P7A28AthCSzVq3g1Vfh7ruh3BmxcNJp2zro27dv3CXL8po3b859993Ht99+i5lx1FFHMW/evJ2OquvqnnvuAWDgwIH89ttvHHXUUbz66qvss88+gHe6e86cOVx33XU89thjdOnShbvvvptLL720dBlHHHEE8+fP55ZbbuGkk06iqKiItm3bcv7554c1VhGJXSWXzM4BngL2DB152mnwz39C4HslkszPtaZ4lJWV5SorQbd06VI6dOhQzxFJrNDnLxK/npk6lXVDh3JFyQ1BQBHw5UUX0emZZ6BB1aUQzGyRcy6rrnHoyFNEROLD999zycMPB++kBX5u0IAlN9/M6X/7W72GouQpIiKx77nnYPBgCK0o1qMH+06Zwr577ln5fBGiG4ZERCR2bdkCV1zhldkrSZypqfDgg/DiixCFxAk68hQRkVj17bdw8cXw2WfBtgMOgJkz4eijoxcXOvIUEZEoqrRa0LPPQufOZRPnhRfCp59GPXGCkqeIiERJRR1sjBw6lG9PPRUuuQRKOrRo1Agefhhmz4ZAr1LRptO2IiISFeWrBR0CzCospN3bbwcnOvBAmDXLOwqNITryFBGRqAitFnQpkAuU6TSsVy/vNG2MJU5Q8hQRkShJT0+nKfAk8E+gWaB9K8CkSV5R9xYtohVelZQ8ZSe33347hx12WJ2Xc/LJJ3PllVeGISIRSUQPXX45uWaEdmf/jRlvjBkDl10GFfTEFCuUPOPQ4sWLadCgQZn+MaNp6tSpNGvWbKf2OXPmMHbs2ChEJCIxb+pUzrn9djqGlIh9YZddWDx5MmfdfHMUA/NHyTMOPf7441x++eUsWbKEpUuXRjucSu2xxx40b9482mGISCzZtAn694eBA6Gw0Gtr2hSefJLzf/+dXkOGRDc+n5Q840xhYSFPP/00Q4cOpWfPnjz55JOl4/Ly8jAznn/+ebp160ZaWhodO3bk9ddfL52mqKiIwYMHc8ABB9C0aVPatWvHvffeW2ln2fPnzyc1NZVVof3k4d0ld8QRR/DOO+8wcOBANm/ejJlhZtx+++3Azqdtt23bxs0330xGRgaNGzembdu2TJgwIYxbR0Ri2pdfes9oTp8ebOvQAT7+GAYNiunTtOUpecaZ5557joyMDI444gj69evH9OnT2b59e5lpsrOzGTlyJJ9//jlHH300vXv3ZlPgeani4mL2228/Zs2axdKlSxkzZgx33XUXTz31VIXrO/HEEznwwAOZHrKzFxcXM336dAYPHsxxxx3HP/7xD9LS0li5ciUrV67khhtuqHBZ/fv3Z/r06YwbN46lS5fy5JNPsluMPLMlIhHkHDz+OBxzDPzvf8H2/v3hk08gDPdY1DvnXEK+unTp4irz1Vdf7dzofbzRedXAiSee6O677z7nnHPFxcUuIyPDPffcc84553744QcHuEmTJpVOv2LFCge49957r9Jljh492p122mmlw7fddps79NBDS4fvu+8+1759+9LhV155xTVq1MitXbvWOefcU0895XbZZZedlnvSSSe5K664wjnn3DfffOMAN2/evBr9vZFQ4ecvInU2Y8YMl5GR4czMZWRkuBkzZji3caNzl1xS9jsvLc25qVOjEiOQ68KQY3TkGUe+++47PvjgA/r06QOAmdG3b1+eeOKJMtMdcUTwSal9990XgDVr1pS2TZo0iaysLFq1akWzZs148MEHyzxvVV7//v35/vvvWbBgAQBTpkzhvPPOY88aFGRevHgxKSkpnHLKKb7nEZH4UVG1oIeHDGHjwQd7j5yUOOww72izf//oBRsG9VZhyMxuAi7AKyKxFfgQuMk5t6SKeTKBHyoYdYZz7tUIhBnTnnjiCYqKikhPTy9tc4E71X788cfSttTU1NL3FriGUHJNc+bMmVxzzTXcf//9HHfccbRo0YKJEyfywgsvVLreVq1a0aNHD6ZMmcIhhxzCSy+9xL///e8axV4Sp4gkpvLVgoYDD27ZQpPQ+yWGDIHx4yEtrf4DDLP6LM93MvAI8AlgwN+AN8yso3NuXTXzdgc+Dxmubvqai/Ev9x07djBt2jTGjh3L2WefXWZcv379eOqpp/jLX/5S7XLef/99/vCHP5S5kWfZsmXVzldyg1Lbtm3Ze++9+dOf/lQ6rlGjRhQVFVU5f+fOnSkuLubtt9+me/fu1a5PROJLydmrFsBkoFfoyGbN4LHHIHDWLBHUW/J0zv05dNjM+gEbgOOB6g5jfnXOrapmmoT28ssvs3btWoYOHbrT6dLevXvz6KOPcumll1a7nIMPPpipU6cyb948DjroIJ599lneffdddt999yrn69atG3vuuSd33HEHN954IykpwTP+mZmZbNmyhddff52jjjqKtLQ00sr9smzXrh0XX3wxQ4YMYfz48XTu3JkVK1aQl5dHv379arAlRCQWpaens2d+PrOAA0Pav0pNpeOiRXDwwdEKLSKiec2zeWD9631MO8fM1pjZB2bWM8JxxaQnn3ySU045pcLrjBdddBH5+fm88cYb1S7nsssu4+KLL6ZPnz4cffTR5OXlcf3111c7n5kxcOBAtm/fzsCBA8uMO+644xg+fDiXXHIJrVq14t57761wGdOnT6dPnz6MHDmS9u3bM2DAADZs2FDtukUkxjnH7JNOYgFlE+fjDRvyxWOPJVziBLBoXYsys1lAOyDLOVfhOT8zawn0Bz4AdgA9gGygv3NuRgXTDwOGAaSnp3fJz8+vcN1Lly6lQ4cO4fgzksqIESP47rvvyjw3Go/0+YuE0W+/weDBMGdOadNG4KaWLTnuH/+gb9++0YutAma2yDmXVdflRKVLMjMbB/wR+GNliRPAObcWeCCkKTeQUEcBOyVP59xkvNPtZGVlxfZFzDiyYcMGFi1axPTp05k1a1a0wxGRWPHJJ3DxxZCXF2zr3JkWM2cy8aCDohZWfaj307Zm9iBwCXCqc+77WiziI7wjVqkn5557LmeffTaDBg3irLPOinY4IhJtzsGECXD88WUT55VXwoIFkOCJE+r5yNPMxgO9gZOdc/+rbvpKdAJWhi8qqc4777wT7RBEJFZs2OCdpn3++WDbrrvCk0/ChRdGL656Vm9HnmY2ERiId9S53sxaB17NQqYZa2Zvhgz3N7M+ZtbBzA4xsxuAK4CH6ituEZFklJOTQ2ZmJikpKWRmZpKTk+N1TN2lS9nE2aWL155EiRPq98jz8sC/b5ZrvwO4PfB+H8rerAVwC5ABFAHfAIMqulmoppxzpQUEJHmoWINI9UqqBZUUPcjPz+ejgQPp7RwNduwITnjFFfDAA9C4cZQijZ76fM6z2kzlnBtQbngaMC3csaSmplJYWLjTs4iS+AoLC8tUYBKRnYVWC2oGPAb0Ce2AonlzeOIJ72ahJJWUtW332msvfvrpJwoKCnQkkiSccxQUFPDTTz+x1157RTsckZhWUi3oMCAXKFMX6MgjYdGipE6cEKVHVaKtRYsWAPz88887decliSs1NZW999679PMXkYqlt2nDqcuXMxFoGtL+dLNm9Fm40Ou8OsklZfIEL4HqS1REpJzNm3krPZ22IT0tbQZGNmrEqZMmKXEGJOVpWxERqcBXX8Exx9D2/fdLm5YAPfbdl1OnTIm5akHRVO2Rp5l1wHu85CQgE+8o/hfgU2Ae8LxzbmsEYxQRkUibMQMuuwxCuhWjf38OmziRN3fZJXpxxahKjzzNrLOZvQEsxuv5ZAFwP3Az3h2wDhgD/Gxmo80s+e5VFhGJd4WFMHQo9OsXTJxNmsCUKTB1KihxVqiqI88XgHuBi5xzlfZ8YmZdgWuB64G7whueiIhEzDffwEUXwRdfBNsOOQRmz4bDD49eXHGgqmue7ZxzE6tKnADOuYXOuYvxjkpFRCTGVFgtaOZMrzpQaOLs08cr9q7EWa1Kjzydc9uqmtHMUp1z2/1OLyIi9a98taBV+flsHjAAQisFNW4M48fDsGGgymu++Lrb1sxGmtmFIcNPAoVm9rWZHRKx6EREpE5CqwUdgNc58rDQxHnggbBwoXezkBKnb34fVRmJd4ctZnYicDFe0YnPKNvfpoiIxJCSakHn4T0i0SV0ZM+eXrWgo46KQmTxzW/y3A/IC7w/B5jtnJuFV9D92PCHJSIi4XBgmzY8iHcH6G6Btm3AbbvvDrNmed2JSY35TZ4bgVaB990I9oyyHWgS7qBERCQM8vNZmJrKNSFNecBpjRtz8EMP6TRtHfhNnq8BjweudR6EVxwB4FDgh0gEJiIidfDvf8NRR9Fy2bLSpn8BPfbfn+FPPqlqQXXkN3legXeduSXQ0zm3LtDeGXgmEoGJiEgt7NgBo0dDjx6wPvCkYcOGcP/9nFtczBc//qjEGQa+CsM75zYCV1XQflvYIxIRkdpZuRJ694b584Nt++/vXdvs2jV6cSUgv4+qdAx9JMXMupnZDDO7ycwaRC48ERHx5Z13vLtmQxPnGWfA4sVKnBHg97Ttk8BRAGa2P96p8z3wTuf+PTKhiYhItYqL4Z574LTTYPVqr80M7rwT5s6Fli2jG1+C8ps8O+A9IgRwEfCRc+5MoB9ejysiIlIPQkvtHdGmDSuysuDGG70kCtCqFbz2GtxyC6So18lI8dsZdgO8R4MATgNeCbxfBuwd7qBERGRnoaX2OgPPrVjB/itWBCc4/nivZu1++0UtxmTh92fJEmCEmZ2AlzxfDbTvB6yNRGAiIlJWSam9oXh9RB4QOvK66+Dtt5U464nf5DkaGAq8AzzjnPsy0N4D+DgCcYmISDlr8/OZBkwGSjpQ3gBcCPDAA5CaGq3Qko7fR1Xmm1kroEW5LsoeAwoqmU1ERMLl66/JTU2l/fbSzqz4HOgJbM/IiFpYycr31WTnXBHQwMz+YGaNA215zrk1EYtORES85zSzssokzqeArsDPaWmMGTMmaqElK7/PeTY3s9nAGrxT7fsF2ieZ2e2RC09EJIlt2wZXXw29esGmTQDsSE3lr3vuyWAz9srIYPLkyaoYFAV+jzzvAfbFK8dXGNI+Fzg/3EGJiCS9H3+Ek06CCROCbQceSMOPP+a+tWspLi4mLy9PiTNK/CbPHsA1zrnPABfSvhRoG/aoRESS2WuvQefO8OGHwbbzzoPcXOjUKXpxSSm/yXN34NcK2psDReELR0QkiRUVwR13QPfusDbwFGCDBnDffTBnDuy2W9XzS73xmzw/wTv6LFFy9HkZ3jVQERGpodBqQUe1acPPRx0Ft98OLvAVu88+3rObN9ygvjdjjN8KQzcD/zGzQwPzXBd4fwxwYqSCExFJVKHVgv4AzF6xgn1DqwWdcgo8/TS0bh21GKVyvo48nXMLgOOARngl+U4Dfga6Ouc+rWpeERHZWUm1oCuB+UCb0JE33eRd91TijFnmnKt+qjiUlZXlcnNzox2GiEiFWpjxONArpG098Bfg3wn6vRwLzGyRcy6rrsup9LStme3hdyHOuXV1DUREJGksXcqnDRty0I4dpU25eF1WOVULigtVXfNcS9nHUipigWnUIbaIiB+zZ8OgQWUS56PAtUCDtDQmq1pQXKgqeZ5Sb1GIiCS67du9fjfHjStt2tGoEX9t3pzx69aRnp7OmDFjVPQgTlSaPJ1z74ZzRWZ2E3ABcAiwFfgQuMk5t6Sa+Q4HHsa7s3cdXjH6O12iXqwVkcSzapVXYm/+/GDbgQfScM4cHjziCB6MXmRSS35r2xaZ2V4VtO9pZn6LJJwMPIJ31+6pwA7gjaqurZpZC+B1YDVwNDAS+Ctwnc91iohE1wcfeNWCQhPnOed41YKOOCJ6cUmd+H3Os7KncxsD2/wswDn35zILNOuH1xXd8cC/K5mtL5AG9HfOFQJLzKwD3nOm43T0KSIxyzl46CG4/nooub5pBn//u3f6NsV3p1YSg6pMnmZWcoTngOFmtilkdAPgBOB/tVx3c7wj3/VVTNMVeC+QOEv8B7gTyAR+qOW6RUTCKicnh+zsbJYvX84h++/PvP33J3PhwuAEe+7pFT04/fToBSlhU92R51WBfw0YQtk6ttuAPGB4Ldc9HvgMWFjFNK2BFeXaVoeMK5M8zWwYMAwgPT29lmGJiNRMaLWgdsDsH38k88cfgxNkZcFzz4EeQ0kYVSZP59wBAGb2NnCBc66qo0TfzGwc8Efgj4FOtqsMo/zslbTjnJsMTAavSEJd4xQR8aOkWtB5wDSgRejIYcNg/Hho0iQ6wUlE+C3Pd0oYE+eDwCXAqc6576uZfBXeEWaokhuXViMiEgN+ys9nLPACwcS5BRgE8NhjSpwJqKoKQxPwHiXZHHhfKefcSD8rM7PxQG/gZOecn2ulC4F7zKyJc25LoK0bXl3dPD/rFBGJqDVreKdJE47fsqW06QfgQmCdTtMmrKpO2x4OpIa8r4yv06NmNhHoB5wHrDezkiPKTc65TYFpxgLHOOdOC4x7GrgNmGpmfwcOBm4E7tCdtiISdR99BD17lkmc8/AeE9iqakEJraoiCadU9L4OLg/8+2a59juA2wPv9wEODFnvBjPrBkzEK/24HngAGIeISLQ4B5MmwdVXe5WDAGfG+BYtuH7DBtpkZKhaUIJTryoiIjVRUAAjRsD06cG23XaDnBw488zoxSW+RLxXlQpW2AuvH8+9KHejkXOuR10DERGJecuWwQUXwBdfBNs6dYLnn4e2baMXl9Q7v+X57gNm4BUm+A34tdxLRCSxvfKK97xmaOIcMAAWLFDiTEJ+jzz/AlzinHsuksGIiMSc4mIYMwZuu8271gnQqJFXem/oUK/kniQdv8UVU/CqAYmIJLycnBwyMzPZ1YzXmjWDW28NJs42beC997ziB0qcSctv8pwMXBrJQEREYkFJqb0m+fl8BJxeGFJa++STYdEiOOaYaIUnMcLvadvdgD6Bx0a+ALaHjvRbJEFEJNZlZ2fTraCA6ZQts/d48+YMff11aOj7PktJYH73go4ET9u2LzcuMZ91EZHkU1TEkPx8bglpKsDrFePZTZsYqsQpAb72hDAVSRARiV3r10PfvmUS5/fA+Xin2zLUU5OE0M8oEZEvv4Tzz/ee4wz4D14PFuuBtLQ0xqjUnoTwlTzN7KWqxqtIgojErVmzYOBAr3JQwH/POYcRn3/Obz/+SEZ6ukrtyU78HnmWL4SQChwJtAHmhDUiEZH6sGMH3Hwz3HdfsG2XXWDqVA7t2ZPq+kuU5Ob3mufAitrN7AHg97BGJCISaWvXQu/e8GZIPxXt2sELL8Chh0YvLokbfp/zrMxjBHtLERGJfYsXe2X2QhPnWWfBxx8rcYpvdU2eh4QlChGRCCmpFpSSksK1LVuy4w9/gPz84AS33govveT1jCLik98bhiaUb8Lre/MMYEq4gxIRCYeSakHbCgp4ELj615DbN5o3hxkzoIfud5Sa83vD0OHlhouBX4BrUfIUkRiVnZ3NLgUFvAycHNL+XcOGHPTJJ3CITp5J7ahIgogkrD3y83kXyAhpmwMM3LGDDUqcUgd1veYpIhKbnn2WD8xKE2cxkA30BHbPyKh8PhEfKk2eZvaEmVW7h5mnr5n1C29oIiK1UFQEN90El1xC00A3YhuAc4C7gKaqFiRhUNVp2xXAF2b2EfASkAusBLYAu+MVi/8j0AvIAy6LaKQiItXZsAH69IFXXilt2ti6NeenpPDOypWqFiRhU2nydM7dbmYTgaF4iXF8uUl+B94ABjrnXotciCIiPnz9NZx7rvdviTPPpEVODm/pMRQJsypvGHLO/YJ3puMuM9sdSAeaAmuBZc45dUcmItH3yitwySWwcWOw7cYb4e9/hwYNoheXJCzfvao459bjdTAgIhIbnIN77/WucZb8lm/aFKZM8crviUSIuiQTkfhUUACDB8Ozzwbb2rSBF1+Ezp2jF5ckBT2qIiJxo6TUXoYZS3bfvWziPOEEyM1V4pR6oeQpInGhpNReen4+nwCHbdsWHDliBLzxBuy1V9Tik+Si5CkicSE7O5t+BQW8CZSkyG3ATXvsAY88Ao0aRTE6STY1Sv1g0KoAACAASURBVJ5m1tLM/mBmjSMVkIjITrZt46b8fCYBqYGm1cCpwD3rdR+j1D9fydPMmpvZLGANsADYL9A+ycxuj1x4IpL0fvkF/vSnMlVYFgFZwAdAenp6dOKSpOb3yPMevITZGSgMaZ8LnB/uoEREAFiyBI45Bt57r7QpBzgBrwRamkrtSZT4TZ49gGucc58BoYURlgJtwx6ViMjcudC1K+TlecNmLO7dm+z0dLaYkZGRweTJk1VqT6LC73OeuwO/VtDeHCgKXzgikvScg/vvh9Gjg4UPmjWDZ57hqLPPJi+qwYl4/B55foJ39Fmi5OjzMrxroCIidbd1KwwcCKNGBRNnZiYsXAhnnx3V0ERC+T3yvBn4j5kdGpjnusD7Y4ATIxWciCSRNWvg/PNhQcjv8RNOgOefh1atoheXSAV8HXk65xYAXYFGwDLgNOBnoKtz7lO/KzOzE83sJTP7ycycmQ2oZvrMwHTlX939rlNE4sAXX8DRR5dNnIMGeYUPlDglBlWbPM0s1cxmAoXOuf7OucOccx2dc5c6576s4fqaAUuAqyl71251ugP7hLzequF6RSTGlJTaO8+MzZ06wfLl3oiUFHjgAXjiCRU+kJhV7Wlb59x2MzsduKmuK3POvQK8AmBmU2sw66/OuVV1Xb+IxIacnByGDR3KyMJCxgApgeub25s0IfX55+HMM6MboEg1/N4wNAe4IJKBVLd+M1tjZh+YWc8oxiEiYfC3m29mUmEhYwl+CS0Dzth9dyVOiQt+bxhaDtxiZicAucDm0JHOuXHhDixgE3ADXiGRHXh3/M40s/7OuRnlJzazYcAwUNURkZi1ahVTly+na0jTO0BPYN0qnWCS+GDOueonMvuhitHOOVfjQglmtgm40jk3tYbzPQL80Tl3RFXTZWVludzc3JqGJSKR9PnncM458OOPpU2TgSuB7UBGRgZ5JUURRCLAzBY557LquhxfR57OuQPquqIw+ggYGO0gRKSGXn4ZevWCzd6JqyLgWuChwGiV2pN4Eo9dknUCVkY7CBGpgQkToEeP0sRJixa8+9e/8lJGBqZSexKHfB15mtmEqsY750b6XE4z4KDAYAqQbmadgHXOueVmNhY4xjl3WmD6/nhncxYDxcA5wBXAaD/rE5Eo27EDrr0WHn442JaZCS+/zKkdO5J3771RC02kLvzeMHR4ueFUoH1gft9FEvB6EXo7ZPiOwGsaMADvGc4Dy81zC5CBd5bnG2BQRTcLiUiM2bgReveGefOCbcceC//6F+y1V+XzicQBv9c8TynfZmZNgCeB93aeo9LlvANYFeMHlBuehpdYRSSeLF/u1aL9MqSOSq9e8NRT0LRp9OISCZNaX/N0zm0BxgDZ4QtHROJebi784Q9lE+ctt8DTTytxSsLwe9q2Mq3wSu6JiMALL0DfvlAYqL6ZmgqPPw79+0c3LpEw83vD0HXlm/CuT/YlUG5PRJJXzowZ5F91FTf+9lvwdNbuu3vJ9KSTohmaSET4PfK8qtxwMfAL8BQwNqwRiUhceXraNLYOHszNRUWlbd+Z8dXNN9NDiVMSlK8KQ/FIFYZE6sFvv/HePvtwwpYtpU3zgfOB5qoWJDEoXBWGan3DkJkdFLjjVkSSUX4+HHdcmcQ5DegGrAOWl3QxJpKAfCVPM7srULAA87yB98zlSjP7QyQDFJEY9Omn3jObS5eWNmXjPay9LTCszhkkkfk98uwLfB14fwZwJHAsMB24OwJxiUismjcPTjwRAj2gFDVsyIBGjbgrZBLVqZVE5zd57g2sCLw/E5jlnPsYr6bzUZEITERi0BNPeL2ilNSo3W03Grz5Jt2mTCFDdWolifi92/ZXvBJ5K4DTgZtC5q+0YpCIJAjn4NZb4e9/D7ZlZHhHoR060BeULCWp+E2ezwNPm9k3wB7Aq4H2TsB3kQhMRGLEtm0wZAj885/Bts6dYe5c2Gef6MUlEkV+k+d1QD6QDoxyzgXO2bAP8GgkAhORGLBhA1x4Ibz5ZrDtjDNg1ixopuJikrz8FobfATxQQfuDYY9IRGLDihVw5plla9QOGQKPPgoN61rZUyS++X1U5aTQR1LMbICZvW9mjwX66BSRBJGTk0P3ffdlRZs2ZRPn3/8OkycrcYrg/27bfwCtAczsEOAx4AugK3BfZEITkfqWk5NDzuDBzFq5kv0DbduBBcOHQ3Y2mO4PFAH/yfNAoOQn6IXA6865y4GhwDmRCExE6t+Ca67hxa1baREY3gB0B/qEdmgtIr5vGHJAg8D704AXAu9XAXuGOygRiYKHHuKhtWtLf1GvwKuIsgQwldoTKcPvkecnwP+ZWT/gBKDkZ2gmXgIVkXjlnHdKduTI0i+Er4Dj8BInqNSeSHl+jzyvAZ4GzgXGOOeWBdovAhZEIjARqQc7dsBll8GUKaVNH6WkcGZxMesCwyq1J7IzX0eezrklzrkjnHO7OufuCBl1A14taBGJNwUFcMEFZRInZ53FD5Mn01yl9kSqVKN7zs0sC+/mobmBQgkNgB2RCExEImjdOujRAz74INjWvz88/ji9U1PpPXhw9GITiQN+n/Pc28w+Aj7GO327d2DUOCooniAiMWzFCjjhhLKJ88Yb4amnIDU1enGJxBG/Nww9SPDO2oKQ9tl4heJFJB4sXQrHHQdffRVse/BBGDtWz3CK1IDf07anAac559Zb2f9gy/Dq3YpIrPvkE68u7a+/esOpqTBtGlxySXTjEolDfpNnU4IdxIdqBWwJXzgiEhHvvOP1w7lpkzfcrBnMmQPdukU1LJF45fe07XzK3lXrzKwBMBp4s8I5RCQ2/Pvf0L17MHHusQe89ZYSp0gd+D3yHAW8a2ZHA43xbhI6FNgVOD5CsYlIHX1w+eX84dFHS/+jF+y+O2nz58Ohh0Y1LpF45/c5z6+Aw/EKIrwGNMG7WeiokIIJIhJDPh4wgK4hiXMZkFVYSM5nn0UzLJGEYM65qicwSwXeB/7inPu6XqIKg6ysLJebmxvtMETqn3Pe3bPZ2aVNX+LdFr8KyMjIIC8vL0rBiUSXmS1yzmXVdTnVnrZ1zm03swPwisOLSCxzDkaNgvvvL236EDgTWB8YXq4i7yJ15veGoWl43Y+JSKwqKvLq1IYkzjeAPxFMnKAi7yLh4PeGoV2AvmbWDVgEbA4d6ZwbGe7ARKQGduyAAQMgJ6e0aXlWFhf9979sLiwsbVORd5Hw8Hvk2QH4FO8HbFu8m4dKXodFJjQR8WX7dujTp0zipH9/0hcu5OHHHydDRd5Fwq7aG4bilW4YkqSwdSv06gX/+lewbfhwmDgRUvz+NhZJHuG6YahG/7vMrKmZHRZ4Na3pyszsRDN7ycx+MjNnZgN8zHO4mb1rZoWB+W41UxFOEbZs8boUC02cV18NjzyixCkSYX57VWlsZv8A1gGfA18A68xsvJk1qcH6muF1Tn81UFjNtJhZC+B1YDVwNDAS+CtwXQ3WKZJ4Cgq8cnuvvBJsGzXKK/Ku35YiEef3hqFH8R4TGwIsDLR1BcYCzYFBfhbinHsFeAXAzKb6mKUvkAb0d84VAkvMrANwnZmNc4l6zlmkKps2wdlnw7vvBtv+7//gjjuUOEXqid/keRFwgXPu9ZC2781sDfA8PpNnLXQF3gskzhL/Ae4EMoEfIrRekdi0YQOceSYsWBBsu/NOuOWW6MUkkoT8XhjZDPxUQftP+Dj9Wget8U7ZhlodMq4MMxtmZrlmlvvLL79EMCyRKFi/3ivmHpo477tPiVMkCvwmz4eA20JvEgq8/7/AuEgqf2rWKmnHOTfZOZflnMtq1apVhMMSqUe//cavXbp4fXIG5PbrBzfcEMWgRJKX3+R5LHAG8JOZvWNm7wAr8Kp+HRO4g/YlM3spzPGtYucjzL0C/5Y/IhVJTBs2sDYriz1/CF6luAw46fnnyQl9tlNE6o3fa55r8a5thqqP640LgXvMrIlzrqTT7W7Az0BePaxfJLp+/x3OOIOWy4KdFw0DHgcoKCA7O1tFD0SiwFfydM4NDMfKzKwZcFBgMAVIN7NOwDrn3HIzGwsc45w7LTDN08BtwFQz+ztwMHAjcIfutJWEt2mTd3PQwoWlTSMIJM4AFXkXiY76fpI6C1gceDUF7gi8/1tg/D7AgSUTO+c24B1p7gvkAhPxOuIeV38hi0TB5s1w1lnw/vulTVcBk8pNpiLvItHh97RtWDjn3iF4w09F4wdU0PYlcGLkohKJMSUFEObPL21a1LcvU154wRsXoCLvItGjGl4isaSwEM49F95+O9h23310mTGDyZMnq8i7SIxQYXiRWLFtG5x3HsybF2wbOxZuvDF6MYkkmKgUhheRCCkqgksvLZs477xTiVMkRvlOnmZ2hpnNNbOvzKxNoG2ImZ1W3bwiUgXn4LLLYPbsYNstt6hykEgM89urSl9gFvAtcACQGhjVABgVmdBEkoBzXpWgJ58Mto0cCX/7W+XziEjU+T3yHAUMdc5dC+wIaf8Q6BT2qESSxZgxMC7kyav+/dWtmEgc8Js82xHsiizUJqBF+MIRSSITJnhdiZU4/3x44gl1ZC0SB/z+L/0Zr7pPeScCyypoF5GqTJsGV18dHO7WDZ55BhrW66PXIlJLfpPnZGCCmR0fGG5jZv2Be/E6yhYRv+bOpXhQsAvc3MaNebZXL2jcOIpBiUhN+K1te6+Z7Qq8DjQB3ga2Avc75yZGMD6RxPLRR+y48EIaFhcD8Bnwp61b2T5yJEVNmqjogUicqLZIgpk1BE4HPsLr+Loj3hHrV865TRGPsJZUJEFizjffwHHHwa+/AvA9cBzBvvUyMjLIy8uLUnAiySFcRRKqPfJ0zu0wszlAe+fcr3gF2kWkJlatgu7dSxPnL0B3ynZKqx5SROKH32uenxPsSkxEauL3370eUgKdWReacTbeQ9Oh1EOKSPzwmzxvBx4ws/PMrI2Z7RH6imB8IvFt+3bo2RM+/dQbbtCAj667jiVpaWUmUw8pIvHFb/J8GTgcmAPk4Z11+gVYG/hXRMpzDkaMgNdeC7ZNmsTJ99+vHlJE4pzfh8pOiWgUIolo3LiyZfduuw2GDAGgb9++SpYicczvoyrvRjoQkYTyr3/BX/8aHO7Xz0ueIpIQfJczMbO9gSvwHlVxwH+BR51zq6ucUSTZLF4Mffp4p20Bjj8eHn9c9WpFEojfXlWOB74D+uA967kFuBT41sy6Ri48kTjz889wzjlQUOANt20LL7yg6kEiCcbvkef9wDPAcOdcMYCZpQCTgAfwnvUWSW4FBXDuufDTT97wrrvC3LnQqlV04xKRsPObPDsBA0oSJ4BzrtjMxgGLIxKZSDxxDgYPhpKqVg0awKxZ0KFDdOMSkYjw+6jKBrxOsMs7APgtfOGIxKlx4+DZZ4PDDz0Ep58evXhEJKL8Hnk+CzxpZqOABXg3DP0RuBvvdK5I8nrrLRg1Kjg8fLj3fKeIJCy/yXMUYMCUkHm243VHdmME4hKJD/n5cPHFEOglha5dYfz46MYkIhHn67Stc26bc+5qYHe8659HAXs45651zm2LZIAiMauwEC64oLTY+0rgmOXLyZk9O7pxiUjEVZk8zSzNzCaa2U9mtgZ4AljpnPvCOVdQPyGKxKgrryytWbsd6Al88tNPDBs2jJycnKiGJiKRVd2R5x3AALzats8C3fBO1Yokt3/+E6ZMKR28Gu9mAICCggKys7OjEpaI1I/qrnleAAx2zj0LYGYzgA/MrIFzriji0YnEoq+/LnND0D/Z+Rel+uYUSWzVHXm2Ad4rGXDOfQzsAPaNZFAiMauw0LtBaPNmAJY1bEhF99Wqb06RxFZd8mwAlL8haAc1qIkrklCuvx6++MJ737gx/7vjDpz65hRJOtUlQQNmmNnWkLYmwONmVnrDkHOuRySCE4kps2fDoyEnaP/xD84aPpzJGRlkZ2ezfPly0tPTGTNmjLobE0lw5kp6fqhopNlTfhbinBsYtojCJCsry+WWlEoTqasff4TDD4cNG7zhiy6CmTPVU4pInDGzRc65rLoup8ojz1hMiiL1rrgYBg4MJs4DDlAXYyJJzm9tW5Hk9fDD8Oab3vuUFO8xlV13jW5MIhJVSp4iVVm6FEaPDg6PGuV1bi0iSa3ek6eZXW5mP5jZFjNbZGYnVDFtppm5Cl7d6zNmSVI7dsBf/gJbtnjDRx4Jd9wR3ZhEJCbUa/I0s17AeOAuvPq4C4B5ZlbdQ3HdgX1CXm9FMk4RAB58MNg/Z6NGMGOG96+IJL36PvK8DpjqnHvcObfUOXcVXj3t6vpv+tU5tyrkpWL0Elnffw+33RYcvuMOOOyw6MUjIjGl3pKnmTUCugCvlRv1GnBcNbPPMbM1ZvaBmfWMSIAiJZyDyy7zqgmBd7r2+uujG5OIxJT6PPJsiVexaHW59tVA60rm2QTcAFwMnAm8Ccw0s0srmtjMhplZrpnl/vLLL+GJWpLP9Onwxhve+5QUeOIJSE2NbkwiElOiUWavfFUGq6DNm9C5tcADIU25ZtYSr3PuGRVMPxmYDF6RhLBEK8ll7Vq47rrg8NVXQ1adn6cWkQRTn0eea4Eidj7K3Iudj0ar8hHQLlxBiZTxf/8H69Z57zMz4c47oxqOiMSmekuegZt8FuH1CRqqG8GuEP3ohHeTkUh4ffEFTJ4cHJ4wAXbZJXrxiEjMqu/TtuOAf5rZx8AHwHC87s0mAZjZWOAY59xpgeH+wHZgMVAMnANcAYzeedEideAcXHONV4oP4PTT4eyzoxuTiMSsek2ezrmZZrYncAve85pLgDOdc/mBSfYBDiw32y1ABt4p32+AQc65na53itTJCy/A22977xs08J7xVO1aEalElb2qxDP1qiK+bdkCHTvCDz94wyNHwvjx0Y1JRCIiXL2qqLatyMSJpYlza7NmHDlnDikpKWRmZpKTkxPl4EQkFkXjURWR2PH773D33aWDN23dyhcrVgCQn5/PsGHDANS5tYiUoSNPSW4TJnjPdgIrGjTg4e3by4wuKCggOzs7GpGJSAxT8pTktX493Hdf6eCtRUVsr2Cy5cuX119MIhIXlDwleU2YABs2eO8PPph30yvu3Ce9knYRSV5KnpKcCgrgoYeCw7feyt/uuou0tLQyk6WlpTFmzJh6Dk5EYp1uGJLkNGUK/Pqr9z4zE3r1om9D779DdnY2y5cvJz09nTFjxuhmIRHZiZ7zlOSzYwe0awd5ed7wQw/BlVdGNSQRqR96zlOktp57Lpg499wTBg6MajgiEn+UPCX5PPJI8P0VV6j4u4jUmJKnJJevv4b33vPeN2gAw4dHNx4RiUtKnpJcnngi+P6cc2CffaIXi4jELSVPSR7btsG0acHhoUOjF4uIxDUlT0keL70Ev/zivd9/f/jzn6Mbj4jELSVPSR6hR52DBnnXPEVEakHJU5LDb7/Bf/4THO7fP3qxiEjcU/KU5PDSS1DSY0qXLtC2bXTjEZG4puQpyWHOnOD7iy6KXhwikhCUPCXxbdsGb7wRHD7//OjFIiIJQclTEt8HH8Dmzd77tm29urYiInWg5CmJb9684Pvu3cEserGISEJQ8pTE99Zbwfd6tlNEwkDJUxJbQQF8/nlw+IQToheLiCQMJU9JbJ9+6vXfCdC+Pey+e3TjEZGEoOQpie3DD4Pvjz02enGISEJR8pTE9tFHwfdKniISJkqekti++ir4vkuX6MUhIglFyVMS26pVwfdt2kQvDhFJKEqekri2b4d167z3KSnQsmV04xGRhKHkKYlrzZrg+5Yt1QWZiISNkqckrtBTtnvvHb04RCThKHlK4lq9Ovi+devoxSEiCUfJUxJXaPLUkaeIhJGSp8SFnJwcMjMzSUlJITMzk5ycnGrnWfzqq6XvJ7/0kq95RET8qPfkaWaXm9kPZrbFzBaZWZXFRs3scDN718wKzewnM7vVLPrdYtTmy7y28yX7PDk5OQwbNoz8/Hycc+Tn5zNs2LAq58vJyWFBSAfY32zcWO08IiK+Oefq7QX0ArYDQ4EOwEPAJiC9kulbAKuAWcBhwIXA78D1PtblMjIy3IwZM1x1ZsyY4TIyMpyZ+ZpnxowZLi0tzQGlr7S0tIjMp3mcy8jIKDN9ySsjI6PKeZ4G5wKvS33MIyKJD8h14chn4ViI75XBR8Dj5dq+BcZWMv0IYCPQNKTtFuAnwKpZV0x9mdd2Ps3jnJlVOI+ZVTnPmyHJs5uPeUQk8YUreZq3rMgzs0ZAAXCJc252SPtE4DDn3EkVzDMd2NM5d1ZI29HAx0Bb59wPVayv9A/LyMggLy+vwukyMzPJz8/fqb2qeVJSUqhou5kZxcXFlYVUq/k0T+0+o8zMTF7Oz+fQwPCRwBfVzCMiic/MFjnnsuq6nIbhCManlkADYHW59tXAnyqZpzWwooLpS8aVSZ5mNgwYVn4h+fn5mNmiStZRYcHTauY5HGhUvtE5t83MvqxkntrOp3lgDyCDstfoi/Pz8/PNbF1l8xxW83kSWUtgbbSDiDJtA20DgEPCsZD6TJ4lyh92WAVt1U1fUTvOucnAZAAzyw3Hr4t4p+2gbQDaBqBtANoG4G2DcCynPu+2XQsU4R0xhtqLnY9GS6yqZHqqmEdERCSi6i15Oue2AYuAbuVGdQMWVDLbQuAEM2tSbvqfgbxwxygiIuJHfT/nOQ4YYGZDzKyDmY0H9gUmAZjZWDN7M2T6p/FuMppqZoeZ2QXAjcA4V/2dTpMjEH880nbQNgBtA9A2AG0DCNM2qLe7bUtXaHY5MArYB1gCXOucmx8YNxU42TmXGTL94cBE4BhgPV6i/ZuP5CkiIhIR9Z48RURE4p1q24qIiNRQ3CbPRKmRWxc12QZmlmlmroJX9/qMOZzM7EQzeynweTozG+BjnoTaD2q6DRJ0P7jJzD4xs41m9ouZ/dvMDvMxX8LsC7XZBom2L5jZFWb2RWAbbDSzhWZ2VjXz1HofiMvkaWa9gPHAXcBReHfrzjOz9EqmbwG8jvd4y9HASOCvwHX1EnAE1HQbhOiOd7255PVWJOOMsGZ4182vBgqrmzgR9wNquA1CJNJ+cDLwCHAccCqwA3jDzPaobIYE3BdOpobbIESi7AsrgNFAZyAL7+940cyOqGjiOu8D4ajxV98v6rFGbqy+arENMvEKS2RFO/YIbY9NwIBqpkm4/aAW2yCh94PA39gM75nyc5J4X/CzDZJhX1gHXBaJfSDujjzNq5HbBXit3KjX8H51VaQr8J5zLvSX+X/wHpPJDHeMkVbLbVBijpmtMbMPzKxnRAKMXQm1H9RRIu8HzfHOqq2vYppE3xf8bIMSCbcvmFkDM+uN9yOisjoCddoH4i55UnWN3PLViEq0rmT6knHxpjbbYBNwA3AxcCbwJjDTzC6NVJAxKNH2g9pIhv1gPPAZXpGVyiT6vuBnGyTcvhC4hrkJ2Ir3WOP5zrnKambXaR+IRm3bcIlYjdw44nsbOOfWAg+ENOWaWUu8Z25nRCa8mJSI+4Fvib4fmNk44I/AH51zRdVMnpD7gt9tkKD7wtdAJ2A3vP6fp5nZyc65JZVMX+t9IB6PPFUjt3bboCIfAe3CFVQcSLT9IFwSYj8wsweBS4BTnXPfVzN5Qu4LNdwGFYnrfcE5t805951zLtc5dxPe0fe1lUxep30g7pKnU43c2m6DinQCVoYrrjiQUPtBGMX9fmBeqc8+eEnjfz5mSbh9oRbboCJxvy+UkwI0rmRc3faBaN8NVcs7qHoB24AhQAe88/ubgIzA+LHAmyHT74r3K+NZ4DDgAry7rK6P9t9Sj9ugP95/rA54/dndEJj/2mj/LXXYBs3w/rN3wquBfGvgfXoS7Qc13QaJuB9MDHyOp+IdSZS8moVMk9D7Qi23QULtC8DdwAl4N/scHvh7i4EzIrEPRP0PrsOGuhzv18FWvKOwE0PGTQXyyk1/ODAf2IL3y+o24vyW9Jpsg8B/lK+AzYEdJBe4NNp/Qx3//pPxrk2Uf01Nlv2gptsgQfeDiv5+B9weMk1C7wu12QaJti8E/r78wPfhGuAN4M+R2gdU21ZERKSG4u6ap4iISLQpeYqIiNSQkqeIiEgNKXmKiIjUkJKniIhIDSl5ioiI1JCSZwIzs5MDndu2DOMyXU17XjCzh83snZDhqWY2N1wxxRozu93MVvvtoDtaavNZxiozGxAoCC61ZGYnmdk3Ztagjsv52MwuCFdcsUrJs54EenZ/o5JxHQJfZOXL7dVk+XlmdkO55gV4ndv+WtvlRsjVgK+eG+LtS9HMDsN70Ho43rafGd2ISpN5ZYWxJYLi7IfifcAYV31B/ercCdxjZgmdXxL6j4sxTwCnmllmBeMG41XGeLOmCw307Vkh5xVJXuVirBKGc26Dc+63aMcRIQcF/n0xsO0Ly09Q1WeWbMwspa5HOsnAzBqamVU/Za2XfxzQHpgdhsW9gtef6BlhWFbsinZJpWR54XX/thK4o1x7Kl4F/1sDwx2Bl4Hf8UpMPQO0Lldiai4wGlgRmOYdypXlCkx7cmC4Zcj8xwJv4ZXk2oCXsPcNjOsOvIfXge46vI5hO5SL1wE9q/g7GwD3B5axHvgH8CjwTvm/IWT4ROBDvNq8G/B6djiMikvP3R6Y51Lgk5DtNBvYL2SZJfOeFlheAV75sc7l4q1qexhe90zLgELgS6ooXwbcXsnnsNNnFmg/HK+EWGFge08Fdq3ks14ViO9uvB+9twf+7lXA6CpiGlDBNhwQ8lkOC2y7zcD35f8+YD+82p8ln+fLQLtq9vUWgc98JV7Zs6VAr5B4NuH1H7kE2BH4rI/G68x9LV6puPeBrjVdbsi0uwMf4O3Du1QSpwHXA9/ilXRbAYz1+7cHPoMlQO/APvI78CKB/28V7Q/AyTVc9oDAsosIqVNbxbb/C96Zpsbl2nOAl6qY72FgTgX78xK8Mn55gc/tKaARXmnQHwPrTE8EaQAACPxJREFUGgeklJt3CjCjvr5fo/GKegDJ9ML74ssP3dHwihEXAW3wTvOtBe7BK9Z8BPBv4OOSefC+UH8P/Gc4DO8LeI/AjnwHgYLQgWlPJiR5AkfifVFPxise3gG4jGAR8QsDr3aBdc8CvgMahcRbXfIchfclfzHeL9mH8L4M3wmZZiqB5In3o2I9XsI9MDBPSbHqRnineDdTrtA1MAjvC7gtcAzwNjA/ZB0lf/vHwCmB5f4H7wvXfG6PMXj9A3YHDgjEtRk4q5K/vRleoX5X7nOo6DNLA37C+7I9HDgJ+AZ4vtx22ojXqW97vK6mioFX8YpcH4x3etgBXSqJqWlg2/4vZBs2DfksV+D9EDkosMxtBDsXSAvENDWwP7THO4OSD6RVsj7DS1pfBbZbW7wjkPMD4wfgJcwFwPGBv6E5XkHzfoHPoD3el/l6gvuun+VuCrzfB++HzixC9t0KYh0L/Ia3Lx0EdAUu9/u34yWXTcALgWm6BsY/FrI/zAReD9n2jWqw7M14Pyg64+03DX18xzQNbLeLQ9p2xfvxeG4V830O3FyureTvmxNY/58Dw/PwkmgH4HxgO3BhuXmHAz9G+zs3kq+oB5BML/6/vXONsauq4vjvP0WtYjWWsQ1SaRAUWjFY0JoALRhbTDCSfiCKSRPbwQIT0URCSJTUV4LSAAlN8AO2KkN4RKCJ1dTSaG2/UPmAI5BSI88pkpJCg6lFytiG5Ye1b+fcM+d5Z+xjXL/k5s7dd5+1n7PXXmuvc48rJQMuy6RtAjanv39M5lf/U9qH0jUL0+d7gNcZv7McAW7MpV1Kt/K8H3i8RX1PxhX7xZm0OuW5B7g587kvLRTbM2n3MKY8ZyaZl5TIW0HGoqgo95wkZ06u7dkfhr4ol6e0P1LbDwKLcul3Ar+vqMeVJIsz196uMQNW4ZuMGQXjdVbmun8A0zJ5ngCerhv73Pc/BHYWpBvdltZJ+CK7PH0ewK0yZfJMw62Nr5SUtRRX8PNKvl9BhbLP5BNuYS5vIfdNXAm+iFuofRXy349br9eVfF/b9tSvb9PtLbgZeL5orvcg+xAwu27uF9T9LuDRzOdB3ENRqnzxTcTKgnlzMNe+R9Jczm6otwN35a69Io1XrcI/UV9x5nkUMbPn8F/wHwCQ9BF8N7c+ZbkAWCzpzc4LXzzBrbIOO81stIcqLKDiXFXSmZIekPSCpH/h7uQ+4PQmwiV9EN/1/7mTZmbv4G7TQsys467cImmTpBskfbRBWedL2ihpt6QDuFKhoK5PZ/7ek947D7yt6o/5wHTg0dx4DNI9Fk3Jj9k8XAkeyKTtwBec+Zm0XdYdwLEXt6rIpc2iN470j5kdxhfGjqwLcIv7QKb9+/ENXVkfLABeNbO/VZR5GH9I8REkzZJ0d4r23I9b6rMYG88mct+Nu3s3m9lgmntlzMef81g2/k3bvtvM9mc+76F+LJrKfsXMenkw9zpgqaQ56fMAMJTGt4z34huBPC/n2rcXeNb8mcLZtHybD+IboOlMUU461hX4P2Q9sE7STHy3/Abw2/RdH26J5qNmofvJ5v/usey6gIPf4a7Ea9P7YdxN9j8NcDGzlZLuxN1xVwC3SFpmZluK8ks6GXfB/hF39b0G9OPntfm6HsoWld47m8aq/ujk+TLwcoXMpuTHTJn65Mmm58uykrReN8JVsvpwJXdVwXVvlMhrEtQyauMjOoeA2cB3GHvM3lbGxrOJ3EO4m/NySXPNbHdF3jp5Tdvey1g0ld3T/7mZPSVpGFgh6TfAZ6iPbt+HK+88TedfPuhrJvC2mZ0wkfJtCeV59HkEPwdcju8I7zWzzmQcxs8Kd2fSmvIfxk/gPMP42dI4JJ2CW0PfNLNtKe18WswRM9sv6VXGgnBIEYILqXk6vZk9hZ+7rJG0GQ9S2EJxu87BleX3zOylVE4v95WV9ge+aRjFz//+1IPsOnYBA5JmZKzPC/GFtcq66oUmc6OIYfycdZ81j44eBk6VNK/GSsxzMfBtM9sEIGk27sVoI9fwDekQsE3SpWaW3/h06IzvF3AXalE72ra9iKK+nyzZVazD4w/6gcfM7O81+f9Kt8djopyLt3PKEm7bo4z5rQsP4OcJZwK/yHz9M/xw/9eSPifpY5KWSPq5pBk1okeARZJOq/hRhNuABUneeZLOlvQNSafjQQb7gFWSzpJ0CR6oUuXqKWItcJOkKyWdjZ8RnlqWWdIZkm6VdKGkuZI+jwdR7Mq0a7qkpZL6Jb0PtwRHgetTH30Jv7esLaX9kRTa7cDtkgZSn3xa0nWSrumhrDz345bFvZI+JWkxcDce8fj8JMjPMgLMTa7ufknvaVHHvcDGdAP9GZIWS7pD0sdLrtmKu+k3SPpiumappGU1ZT0LLJc0X9Jn8UjUrGuwkdzkqv067gLfnub2ONL4rgV+KmllOrJYKGlwAm0vYgQ4N82tfknvmkTZVTyIBygN0r3GlLEF38BMFovwwLYpSyjPY8N63EWyI7uLNrM9eFBLJ6LyGVyhjqZXFd/HI3ZfwM+txmFmTwJLcMvtcXwxugo4lBadr+KKa2cqd3WDcvPcgUfirU/y+/DFooy38IjLh/EFdCjlX5PqvANX4g+mdt1kZq/jC+QyXMn+ALihZT0r+yNlWY1vcm7Ex+IPeDTyS23LKij7Lfy8+wN4RPBG/Kx4YKKyC9iA33u3Fe/Dr7Wo42I8AOdhPGJ3CJ+7/yy55h08CvYx4D7cil5Lvet/AA/i+QuuOH+JK57WcnMKdFuZAgW+i8+z1UneBmBOr20vYV2S/QTe9xdNRLb8R0NMxfeLHyFtDh7CNyAPNajnfcAnJH2yQd5KJJ2Ge1F+NVFZxzOdkP0gCILgOEfSj/CI7vNqAoBIxx+vmNmqhrJvBT5sZldPsI634RG6k+GhOW6JM88gCIITh8uB66sUZwpGXAJcht/L3JSfAN+SNK0goKsNr+FHHlOasDyDIAimEJJG8GjXW8xszTGuzpQllGcQBEEQtCQChoIgCIKgJaE8gyAIgqAloTyDIAiCoCWhPIMgCIKgJaE8gyAIgqAloTyDIAiCoCX/BVTcTQAgWIlwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_TFVM = []\n",
    "y_TFVM = []\n",
    "\n",
    "for i in range(0, n):\n",
    "\t\n",
    "\tp_TFVM.append(p[int(n/2)][i])\n",
    "\n",
    "\ty_TFVM.append(y[int(n/2)][i])\n",
    "\n",
    "\n",
    "'''\n",
    "y = []\n",
    "p = []\n",
    "for line in open('pressure.curve', 'r'):\n",
    "\tif not (line.strip().startswith(\"#\") or line.strip()==''):\n",
    "\t\tvalues = [float(s) for s in line.split()]\n",
    "\t\ty.append(values[0]-5.)\n",
    "\t\tp.append(values[1])\n",
    "'''\n",
    "\n",
    "fig = plt.figure(figsize=[7.,5])\n",
    "\n",
    "plt.subplot(111)\n",
    "#plt.plot(y,p,'ko',label='GEOSX')\n",
    "plt.plot(y_TFVM,p_TFVM,'ko',label='This model')\n",
    "\n",
    "r = np.linspace(1,5,500)\n",
    "plt.plot(r,analytic(1, 1, r, np.pi/2.0),'r',linewidth=3,label='Analytic')\n",
    "\n",
    "plt.xlabel('Vertical distance from the crack center, y (m)')\n",
    "#plt.xlim(-3.,3.)\n",
    "plt.ylabel('Pore pressure (unitless)')\n",
    "plt.xlim(0,3)\n",
    "plt.ylim(0,3)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.008526,
     "end_time": "2023-10-06T07:42:43.757460",
     "exception": false,
     "start_time": "2023-10-06T07:42:43.748934",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "papermill": {
   "duration": 18.681246,
   "end_time": "2023-10-06T07:42:43.877079",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2023-10-06T07:42:25.195833",
   "version": "2.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}