MPM versus FEM simulation of beam bending

FEM_beam_bending.edp

 // Macro
 //Macros for the gradient of a vector field (u1, u2)
 macro grad11(u1, u2) (dx(u1)) //
 macro grad21(u1, u2) (dy(u1)) //
 macro grad12(u1, u2) (dx(u2)) //
 macro grad22(u1, u2) (dy(u2)) //
  
 //Macros for the deformation gradient
 macro F11(u1, u2) (1.0 + grad11(u1, u2)) //
 macro F12(u1, u2) (0.0 + grad12(u1, u2)) //
 macro F21(u1, u2) (0.0 + grad21(u1, u2)) //
 macro F22(u1, u2) (1.0 + grad22(u1, u2)) //
  
 //Macros for the incremental deformation gradient
 macro dF11(varu1, varu2) (grad11(varu1, varu2)) //
 macro dF12(varu1, varu2) (grad12(varu1, varu2)) //
 macro dF21(varu1, varu2) (grad21(varu1, varu2)) //
 macro dF22(varu1, varu2) (grad22(varu1, varu2)) //
 
 //Macro for the determinant of the deformation gradient
 macro J(u1, u2) (
       F11(u1, u2)*F22(u1, u2)
     - F12(u1, u2)*F21(u1, u2)
 ) //
 
 //Macros for the inverse of the deformation gradient
 macro Finv11 (u1, u2) (
       F22(u1, u2) / J(u1, u2)
 ) //
 macro Finv22 (u1, u2) (
       F11(u1, u2) / J(u1, u2)
 ) //
 macro Finv12 (u1, u2) (
     - F12(u1, u2) / J(u1, u2)
 ) //
 macro Finv21 (u1, u2) (
     - F21(u1, u2) / J(u1, u2)
 ) //
 
 //Macros for the square of the inverse of the deformation gradient
 macro FFinv11 (u1, u2) (
      Finv11(u1, u2)^2
     + Finv12(u1, u2)*Finv21(u1, u2)
 ) //
 
 macro FFinv12 (u1, u2) (
       Finv12(u1, u2)*(Finv11(u1, u2)
     + Finv22(u1, u2))
 ) //
 
 macro FFinv21 (u1, u2) (
       Finv21(u1, u2)*(Finv11(u1, u2)
     + Finv22(u1, u2))
 ) //
 
 macro FFinv22 (u1, u2) (
       Finv12(u1, u2)*Finv21(u1, u2)
     + Finv22(u1, u2)^2
 ) //
 
 //Macros for the inverse of the transpose of the deformation gradient
 macro FinvT11(u1, u2) (Finv11(u1, u2)) //
 macro FinvT12(u1, u2) (Finv21(u1, u2)) //
 macro FinvT21(u1, u2) (Finv12(u1, u2)) //
 macro FinvT22(u1, u2) (Finv22(u1, u2)) //
 
 //The left Cauchy-Green strain tensor
 macro B11(u1, u2) (
       F11(u1, u2)^2 + F12(u1, u2)^2
 )//
 
 macro B12(u1, u2) (
       F11(u1, u2)*F21(u1, u2)
     + F12(u1, u2)*F22(u1, u2)
 )//
 
 macro B21(u1, u2) (
       F11(u1, u2)*F21(u1, u2)
     + F12(u1, u2)*F22(u1, u2)
 )//
 
 macro B22(u1, u2)(
       F21(u1, u2)^2 + F22(u1, u2)^2
 )//
 
 //The macros for the auxiliary tensors (D0, D1, D2, ...): Begin
 //The tensor quantity D0 = F{n} (delta F{n+1})^T
 macro d0Aux11 (u1, u2, varu1, varu2) (
       dF11(varu1, varu2) * F11(u1, u2)
     + dF12(varu1, varu2) * F12(u1, u2)
 )//
 
 macro d0Aux12 (u1, u2, varu1, varu2) (
       dF21(varu1, varu2) * F11(u1, u2)
     + dF22(varu1, varu2) * F12(u1, u2)
 )//
 
 macro d0Aux21 (u1, u2, varu1, varu2) (
       dF11(varu1, varu2) * F21(u1, u2)
     + dF12(varu1, varu2) * F22(u1, u2)
)//
 
macro d0Aux22 (u1, u2, varu1, varu2) (
      dF21(varu1, varu2) * F21(u1, u2)
    + dF22(varu1, varu2) * F22(u1, u2)
)//

///The tensor quantity D1 = D0 + D0^T
macro d1Aux11 (u1, u2, varu1, varu2) (
      2.0 * d0Aux11 (u1, u2, varu1, varu2)
)//

macro d1Aux12 (u1, u2, varu1, varu2) (
      d0Aux12 (u1, u2, varu1, varu2)
    + d0Aux21 (u1, u2, varu1, varu2)
)//

macro d1Aux21 (u1, u2, varu1, varu2) (
      d1Aux12 (u1, u2, varu1, varu2)
)//

macro d1Aux22 (u1, u2, varu1, varu2)(
      2.0 * d0Aux22 (u1, u2, varu1, varu2)
)//

///The tensor quantity D2 = F^{-T}_{n} dF_{n+1}
macro d2Aux11 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT11(u1, u2)
    + dF21(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux12 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT11(u1, u2)
    + dF22(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux21 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT21(u1, u2)
    + dF21(varu1, varu2) * FinvT22(u1, u2)
)//

macro d2Aux22 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT21(u1, u2)
    + dF22(varu1, varu2) * FinvT22(u1, u2)
)//
 
///The tensor quantity D3 = F^{-2}_{n} dF_{n+1}
macro d3Aux11 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//
 
macro d3Aux12 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//

macro d3Aux21 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
     + dF21(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

macro d3Aux22 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
     + dF22(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

///The tensor quantity D4 = (grad w) * Finv
macro d4Aux11 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad11(w1, w2)
    + Finv21(u1, u2)*grad12(w1, w2)
)//

macro d4Aux12 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad11(w1, w2)
    + Finv22(u1, u2)*grad12(w1, w2)
)//

macro d4Aux21 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad21(w1, w2)
    + Finv21(u1, u2)*grad22(w1, w2)
)//

macro d4Aux22 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad21(w1, w2)
    + Finv22(u1, u2)*grad22(w1, w2)
)//
//The macros for the auxiliary tensors (D0, D1, D2, ...): End

//The macros for the various stress measures: BEGIN
//The Kirchhoff stress tensor
macro StressK11(u1, u2) (                   
      mu * (B11(u1, u2) - 1.0) + 2 * la * (J(u1,u2)- 1)*J(u1,u2)
)//

//The Kirchhoff stress tensor
macro StressK12(u1, u2) (
      mu * B12(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK21(u1, u2) (
       mu * B21(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK22(u1, u2) (
      mu * (B22(u1, u2) - 1.0) + 2 * la * (J(u1,u2)- 1)*J(u1,u2)
)//

//The tangent Kirchhoff stress tensor
macro TanK11(u1, u2, varu1, varu2) (
      mu * d1Aux11(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux11(u1,u2,varu1,varu2)
)//

macro TanK12(u1, u2, varu1, varu2) (
      mu * d1Aux12(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux12(u1,u2,varu1,varu2)
)//

macro TanK21(u1, u2, varu1, varu2) (
      mu * d1Aux21(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux21(u1,u2,varu1,varu2)
)//
macro TanK22(u1, u2, varu1, varu2) (
      mu * d1Aux22(u1, u2, varu1, varu2)  + 2 * la * (2*J(u1,u2) - 1) * J(u1,u2) * d2Aux22(u1,u2,varu1,varu2)
)//
//The macros for the stress tensor components: END

// Parameters
//real mu = 5.e2; //Elastic coefficients (kg/cm^2)
real D = 1.e3; //1.e3; //(1 / compressibility)
real Pa = 0.0; //Stress loads   -3.e2

 real E = 500;   // rajouté...
 real nu = 0.35; //0.4 -0.2;   //rajouté...
 real mu = E/(2*(1+nu));
 real la = E * nu / ((1 + nu) * (1 - 2 * nu)) ;
//la=0.0;

cout << "mu= "<<mu << " lambda= "<< la <<endl;
 
real tol = 1.e-4; //Tolerance (L^2)

int m = 40, n = 20;
real f = -15;    

mesh Th = square(80, 40, [0.6*x, 0.1*y + 0.8]);
cout << "premier plot"<<endl;
plot(Th, wait=true, cmm="TOUT PREMIER DESSIN");

// Fespace
fespace Wh(Th, P1dc);  //P1dc
fespace Vh(Th, [P1, P1]);
Vh [w1, w2], [u1n,u2n], [varu1, varu2];
Vh [ehat1x, ehat1y], [ehat2x, ehat2y];
Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress
Vh [ef1, ef2];

fespace Sh(Th, P1);
Sh p, ppp;
Sh StrK11, StrK12, StrK21, StrK22;
Sh u1, u2;

// Problem
varf vfMass1D(p, q) = int2d(Th)(p*q);
matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);
 
p[] = 1; //1
ppp[] = Mass1D * p[];

real DomainMass = ppp[].sum;
cout << "DomainMass = " << DomainMass << endl;

varf vmass ([u1, u2], [v1, v2], solver=CG)
    = int2d(Th)( (u1*v1 + u2*v2) / DomainMass );

matrix Mass = vmass(Vh, Vh);

matrix Id = vmass(Vh, Vh);

//Define the standard Euclidean basis functions
[ehat1x, ehat1y] = [1.0, 0.0];
[ehat2x, ehat2y] = [0.0, 1.0];

real ContParam, dContParam;
 
problem neoHookeanInc ([varu1, varu2], [w1, w2], solver=LU)
    = int2d(Th, qforder=1)(
        - (
              StressK11 (u1n, u2n) * d3Aux11(u1n, u2n, varu1, varu2, w1, w2)
            + StressK12 (u1n, u2n) * d3Aux12(u1n, u2n, varu1, varu2, w1, w2)
            + StressK21 (u1n, u2n) * d3Aux21(u1n, u2n, varu1, varu2, w1, w2)
            + StressK22 (u1n, u2n) * d3Aux22(u1n, u2n, varu1, varu2, w1, w2)
        )
        + TanK11 (u1n, u2n, varu1, varu2) * d4Aux11(w1, w2, u1n, u2n)
        + TanK12 (u1n, u2n, varu1, varu2) * d4Aux12(w1, w2, u1n, u2n)
         + TanK21 (u1n, u2n, varu1, varu2) * d4Aux21(w1, w2, u1n, u2n)
        + TanK22 (u1n, u2n, varu1, varu2) * d4Aux22(w1, w2, u1n, u2n)
    )
    + int2d(Th, qforder=1)(
          StressK11 (u1n, u2n) * d4Aux11(w1, w2, u1n, u2n)
        + StressK12 (u1n, u2n) * d4Aux12(w1, w2, u1n, u2n)
        + StressK21 (u1n, u2n) * d4Aux21(w1, w2, u1n, u2n)
        + StressK22 (u1n, u2n) * d4Aux22(w1, w2, u1n, u2n)
     )
     //Choose one of the following two boundary conditions involving Pa:
     // Load vectors normal to the boundary:
    - int1d(Th, 1)(              //1
          Pa * (w1*N.x + w2*N.y)   
    )
        - int2d(Th)(   //rajouté... gravité
         f*w2
    )
    //Load vectors tangential to the boundary:
    //- int1d(Th, 1)(
     //    Pa * (w1*N.y - w2*N.x)
     //)
     + on(4, varu1=0, varu2=0)  //2
    ;

//Auxiliary variables
matrix auxMat;

// Newton's method
 ContParam = 0.;
dContParam = 0.01;

//Initialization:
[varu1, varu2] = [0., 0.];
[u1n, u2n] = [0., 0.];
real res = 2. * tol;
real eforceres;
int loopcount = 0;
int loopmax = 400;//45;

// Iterations
while (loopcount <= loopmax && res >= tol){
//while (res >= tol){      
    loopcount ++;
    cout << "Loop " << loopcount << endl;

    // Solve
    neoHookeanInc;

    // Update
    u1 = varu1;
    u2 = varu2;

    // Residual
    w1[] = Mass*varu1[];
    res = sqrt(w1[]' * varu1[]); //L^2 norm of [varu1, varu2]
    cout << " L^2 residual = " << res << endl;

    // Newton
     u1n[] += varu1[];

    // Plot
    plot([u1n,u2n], cmm="displacement");
}

// Plot
plot(Th, wait=true, cmm="PREMIER PLOT");

// Movemesh
mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);

// Plot
plot(Th1, wait=true, cmm="DEUXIEME PLOT");
//plot([u1n,u2n], cmm="TROISIEME PLOT");

// write mesh vertices to file 
ofstream ff("fichier_plot_beam_FEM_neohookean.txt");
int NbVertices = Th1.nv;
cout << "Number of vertices = " << NbVertices << endl;
for (int i = 0; i < NbVertices; i++)
      ff << Th1(i).x << " " << Th1(i).y << endl;

MPM-FEM-beam.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "818f56d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "611897a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_MPM = open(\"fichier_plot_beam_MPM.txt\", \"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3cf8ef83",
   "metadata": {},
   "outputs": [],
   "source": [
    "points_MPM_X = []\n",
    "points_MPM_Y = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0ce8e97a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lines_MPM = f_MPM.readlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fb319f3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for line in lines_MPM:\n",
    "    lu = line[:-1]\n",
    "    liste_lue = lu.split(\" \")\n",
    "    points_MPM_X.append(float(liste_lue[0]))\n",
    "    points_MPM_Y.append(float(liste_lue[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "89aa5423",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_FEM = open(\"fichier_plot_beam_FEM_neohookean.txt\", \"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b96d6482",
   "metadata": {},
   "outputs": [],
   "source": [
    "points_FEM_X = []\n",
    "points_FEM_Y = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f09d0ae7",
   "metadata": {},
   "outputs": [],
   "source": [
    "lines_FEM = f_FEM.readlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "06462731",
   "metadata": {},
   "outputs": [],
   "source": [
    "for line in lines_FEM:\n",
    "    lu = line[:-1]\n",
    "    liste_lue = lu.split(\" \")\n",
    "    points_FEM_X.append(float(liste_lue[0]))\n",
    "    points_FEM_Y.append(float(liste_lue[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4b4d0ec5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f305609cf10>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHiCAYAAAA6dsw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABG40lEQVR4nO3df5icd13v/9c7u0nZTWt3YwvHNilFhWJBONJAqkKteM4BSrWCiA2d5ceXzKZFpVweTUElYvH4NTlensoR2sxWbNmlARUOIEWRozRBpYVWaaFisfxMAC00SW26adMk7/PH574z99xzz8w9P++Z2efjuubKzs4999wzk2Re8/68P5/b3F0AAABFWVX0AQAAgJWNMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEXTEzM40s/vM7AnR9dvMbEuDbc81MzezyR487gfN7MXd7mcYmdn5ZnZnD/f3NjNbin5+kpl90cxOabL9TWb2O716/F5JPo8BPNbJ18DMXmBm9w3icYGVjjAyoszsa2Z21MzOSP3+c9EH/7nR9Zui7Q6b2QEz+4SZPT267W3Rtm9M7eNN0e/f1uQQ3izpT9z90d4+s5Z+T9L/6PeDRK/vkeh1iy9nJYLV4dTlF6L73RTd/jOp/V0X/f61TR727ZJ+vx/Px93/XdInJc33Y//jyN0/5e7nDeKxor8bj6T+Tm1rcx+vNbPjqX1cnLh9nZn9n+hxvm5mr0rd/6fM7F/MbNnMPmlmT+7Ns+sNMzvFzN5tZv9hZv9mZr/SZNufNLPPm9khM3swet5nJ25P/r8YXyYG80yQhTAy2r4qaXN8xcx+WNJUxnY73f1USeslPSDppsRtX5L0mtT2r45+nyn6dv0aSQP5tprk7p+R9D1mtnEAD/fT7n5q4vKtxG0zqdven7it5jWNKkI/L+nLjR7IzL5P0k9K+lCD27uuKkl6r6StPdgP+uPZqb9TOzvYx6dT+7gtcds7JR2V9CRJV0i63syeIUnRl5oPSnqrpHWS7pT0fg2Xt0l6qqQnK/xb2dakSvrPkl7k7jOSzpL0r5KuT22zM/VaHe/PYSMPwshoW1QIDrHXSHpPo43dfVnSLZKemfj1ZyVNJ/5TeoZCoPlsk8fdJOmQu+9P/f4HzOwzZvaQmX3YzNZl3TmqOvyXxPWaMryZXWhm/xB9q7k7+e0ucpuklzbY9xfN7NLE9Ukz+66ZPcfMnmBmS9E3pUNm9lkze1KT59mpv5D042Y2G11/saR7JP1bk/v8V0n/mKw0Ra/TNWZ2j6RHoufS8LUxs6eY2R4ze9jMPiHpjNRj3CHp+1t84z0jqp49HO3r5LZm9vTotgMWhuhembjtpWb2T9G31n2WqKolqkmvi247aGZXmtlzzeye6Ln8UZNjkqQnmNn7o+P6RzN7dmL/Z5nZB8zsO2b2VUtU+qK/W39qZu+J7ntvMsia2Y9E+3vYzN4v6QmJ2y42s/2J618zs1+Njvmh6HiS228zs2+b2bfMbEv0nH+wxfMaCDNbK+nnJL3V3Q+7+99J+oikuWiTl0u6193/LPo7+DZJz7aoippj/x69p/8avb/vNDOLbkv/++502PbVkt7u7gfd/YuSFiS9NmtDd//31JeH45KG4r1ANsLIaLtdoUrwQxZKjL+gJtUKMztV4RvRP6VuSoaapoEm8sOSssbSXy3p/1P4JnJM0jtaPYGMYzxb0q2SfkfhG9qvSvqAmZ2Z2OyLkp6dcXdJ2q1EtUjSiyR9193/UeG5nS5pg6TvlXSlpCPtHmMOjyr8R395dP3V6vw13awQvGYUvtE2e21ukXSXQgh5u1IVL3c/Jul+NX7tpPD34+3RPj6nUE2JP8w+ET3GE6PjelcUXiXpkeh5zkTHe5WZ/Wxq35sUvtn+gqTrJP2GpP8i6RmSXmlmP9HkuC6T9GfR875F0ofMbLWZrVIIf3dLOlvST0l6k5m9KHHfn5H0vujYPiLpj6LntEahErUY7ffPFD6wm3mlQrh8iqRnKfowtPAN/Vei5/ODkpo9l7aY2auiwNbock5i8x+JwveXzOytiQ/8p0k67u7JiufdCq+9oj/vjm9w90cUKnnPUH6XSnquwt+vVyr828vz/N7V5LndE20zq/D/yt2JuyaPP2u/55jZIYV/478qKV1pekMUrO8ys1bvO/qMMDL64iDxXyX9i6RvZmzzq9E/yvslnar6bxNLkjab2WqFD9BWwy8zkh7OOhZ3/0L0H9lbFT5g2h2HLUn6mLt/zN1PuPsnFErGlyS2eTg6hiy3SPoZM5uOrr8q+p0kPa4QQn7Q3Y+7+13u/h9NjuVDif8UP5S67bup/zR/KHX7eyS92sxOV/hgSt8/bUbZr+k73H2fux9Rk9cm+kB6rsI338fcfa/Ch3Ras9dOkm51973u/phCWPhRM9ug8EHzNXf/E3c/FoW7D0h6hSS5+23u/vnouO5RCIXpD+S3u/uj7v7XCuFlt7s/4O7flPQpST/S5Ljucvc/d/fHJf2BQgXjwug5n+nu17r7UXf/isI35ssT9/276DU7rvDvJQ5jF0paLek6d3/c3f9czSuCUng/vuXuBxRe3/8c/f6VCj1U90YVyN9usZ8s/5j6O/UiSXL3W9x9psnlG9H99ypUPZ+oEKo2S/q16LZTJT2UeryHJJ2W8/Y8fs/dD0XH80lVX5um3P0NTZ7bsxLHFx9TruNz929EwzRnSPpNhf8fY+9QCMZPVPi/6iYz+/E8x4v+6MU4NIq1qPCf0FPU+Nv377v7bzbagbt/w8zul/S7kv7V3fdFFdZGDir7P4F9iZ+/rvAffXqooJUnS/p5M/vpxO9WK/znFjtN0qGsO7v7/Wb2RUk/bWZ/ofCtOP6QW1SoirzPzGYUQtdvRB9wWX7W3f9vg9vOiCoNmdz976KKxW9K+qi7H+nBa9rstTlL0sEoCMa+rvB8kxq+dunHc/fDZnYg2veTJW2KQm1sUuE1lZltUmgufqakNZJOUag0JP174ucjGddPVWPJ4zoRDZ+cJcklnZU6rgmFcBNLDo8tKwz5TEb3/6bXni30602OIWtfZ0U/n6UQDOuOtw3Pcff7O7ifJCkKYrHPm9m1CmHk/5d0WNL3pO7yPaoG4Fa355F+bZq9n+06HP35PQqVx/jnlsfn7gfM7GZJd5vZ2YkwHfuYmb1XYajq73t4zGgDlZER5+5fV2hkvUShAa1T75H039V6OEEK/Q9Py/h98oPvHIVKxHcztntE0nTi+n9K/LxPocKS/Ha01t1/L7HND6m2XJsWD9VcJumf4//go2+/v+3u50v6MYVv+69uvJuuLan71zT5Qdnstfm2pNloOCWWLN/HTbA/qOav3cn3MBrWWyfpW9Fj70k99qnuflW0+S0KQyAb3P10STdIapq+2pQ8rlUKzdjxcX01dVynufsljXaU8G1JZ1ttSjyn0cY59rU+63i7ZWZXWP3sreSl0TG7qu/BlyRNmtlTE7c/W9K90c/3KjF8F/09+oHE7d1o9u9dZnZDk+d2ryS5+0GF1zg5xJg8/lYmFaog6cAVS75WKABhZDy8XtILU9+K2/V+Sf9N0p/m2PYzkmYsMVUuUrKwVsa0pGsl/blnd6h/TtLl0Zj/RkWl/siSQlXjRWY2YaHp9GIzS/5H/xOS/rLJ8b0vei5XqTpEE0/3++Fo6Og/FMJSPzvo36EwfLY3x7afkPQcSzREZmj42kSh9E5Jv21ma8zs+ZJ+OnX/5ykMtTT79n+JmT0/6qd4u6Q73H2fpI9KepqZzUXv22oLDajx8NRpkg64+6Nm9jyF4bFeusDMXh4FqjdJekyhZ+ozkv7DQqPvVPS6PNPMnptjn59W6G16o4Xm4JcrvEad+FNJr7PQvzUtaXvyRgvTbr/WyY7d/b1eO+sjfflG9Bgvsagh20Lj6VslfTjaxyMKX1auNbO10ZDEZYoqW5L+j6RnmtnPRX8Ht0u6x93/pdvjV/j3flHUw3G6pLeknt+VTZ5bsifkPZJ+08xmo+dXVu3MwJOivyvnmdmqqEL5B5L+KRpek5m9wsxOjW7/bwpDoB/p8PmhBwgjY8Ddv+zuXS2W5e5H3P3/euhNaLXtUYX/BEqpmxaj3/+bwpj+G5XtrQrfug4qjK2fDAzRB99lkn5d0ncUvvn+mqK/q9GHzCMepvg2Or5vK3zQ/Jhqpyf+J0l/rhBEvihpjzqfnnwo9Q2ubs0Ddz/g7n+TGgZodMz/LulvFZ57o22avjYKAWCTpAOSfkv1FZkrFCoWzdwS3feApAui+8jdH1YIeJcrVCT+TdIOheEYSXqDwgfdwwofZHlCbTs+rND4elBhBsjLo0rXcYXQ9Z8VKoTflXSjQqNyU9Hf45cr9FAdjPbfUXXR3f9SIXx+UqE369PRTY9Ff25Q6yGAu1N/p65r8zB+StI9ZvaIpI8pPJffTdz+BoWZcg8oVA+vcve48vAdhT6T/6HwWmxSbd9NnuPP5KG36f0K1b+7FIJtJ35Loan26wr/dv+nu/9VfGP0mr0gunq2pL9SGMb5vKQTkl6W2NfVCv11hyT9T0llr50GjQGzHP9PAnWibxufkvQjeQJMDx/3A5L+2N0/NqjHHBQzO1/SzZKelyfAtLnvJyr8B/4jPviF6lacqGL0BUmnuPsxM/trSVd7mJI6ckb9+DH8CCMA0ANm9jKFqddrFULlCXf/2UIPChgRLYdpLCy/+4CZfaHB7WZm7zCz+y0sBvSc3h8mAAy9rQrDZ19W6EW6qvnmAGItKyNmdpHCtKr3uPszM26/RNIvK8zm2CTpD919Ux+OFQAAjKGWlREPiycdaLLJZQpBxd39doVZFt/XqwMEAADjrRezac5W7QI/+6PfAQAAtNSLFVizForJHPsxs3lFpzBfu3btBU9/eq5zMAEAgCF31113fdfdz2y9Zb1ehJH9ql1tMF4ZsY67VyRVJGnjxo1+551dLY0BAACGhJm1Op1CQ70YpvmIwgnBzMwulPRQtOgUAABASy0rI2a2W9LFks6wcHKq31I4OZfc/QaFlf4uUVh1cFnS6/p1sAAAYPy0DCPuvrnF7S7pF3t2RAAAYEXh3DQAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAAClVcGLnnHmlh4eTVhQVpw4aaXwEAgBXA3L2QBz7XzvQH9DUd0droNy7JEn9Kq1dLp54q7dghlcuFHCYAAMjBzO5y942d3Lewysh+nZ0IIlIcQKp/So8/Lh08KM3PS2auCXtcZi4zaW5uoIcLAAD6ZEh6Rk5IOqFVelyhMqLEnzHTCa1WHFaWlhQFkxMnA8rEhGQmwgoAACOk0DAyq4OqVCSfOlWuCR2fOl3uJnepMvmLmtJhhaAiZYcUU3gKIaCcOFG9JYSV2sCyZo20bh19KQAADJPCwshxTWqtHmnYC1JefZOWdVoIKlE4Wa99qkz+okolKa6mVC8xV6PAUh328TDsMyGtWkUVBQCAIhUWRlbrqLZP/m7u7curb9I+naPy6pu0uCj57BlyTchtUl75Y1UqCmFF8yppMbpXMqykA4rpxAnJXVpa8pNDPWbS2rVUTwAAGJTCwsiz9HmVV9/U+Q4efTT8+YQnSOWyylrQPnuyyrpRizNXy706/ONTp6pSMc1MPapJParacBLP3qk2zi4vh+rJmjXS9DRDOwAA9NOQNLCmLCxUw8Ypp+S7zzXXhDKHmbRzZ93N5bJ08Aln6XFNyW1SldKntH69tEm3qz6cSPGwzpEjyaEdqiYAAPTacIaRFsEiU6pS0mqb8uJF2rd9Qbfbj4fqycw6eWVBFc1HjbP1s3mkatVkwh7XKnP6TQAA6NJwhpE8waKVPNWVdOi55hqVdaOWdZoqU2/SzIw0pcOJoR0pHtY5odVy2ckpxoQTAAA6M5xhJK0XwzZZ+0iHnvi6pPIfPlMHdy5oWaeFoZ2pU7V+vVRdJVZKDumkw8kae0zr1j7KkA4AAC0MRxiJw0HeCkYe6aCRFU6OHMm+79RU9T5SuM8f/qG2X7xX67VPJS1qvfZrejpues3oN9EpOrj8BM3PS3MXfokT7wAA0MBwhJFLLw3Lp156afbt6WDRSaUkK5zETjklV+WkfOvPap/O0aK9Vvsqf6nrXl4bTjZtMtVPI5aW7niqbP/XZfOvD8M4nBUQAICTCjtR3kYzv3NmJkxVWbcu/Dk7Kx04EObTHjkSKhTLy/XX4+3NpF27wg63bg2Vj3ifrfYRX5ekSiWEk3b3GR+HVLPNwpErtE07dHhynY4dq574L6hWUUpTH9Ti8iv6/EoDANB/I3miPEnVSkhcgUj0bDSVZwgmb+UkHpLpZJ/J401sU9aNOjjzA3r8calk71XtKrHxmiartHTk58J5dBjGAQCsYMVWRvJWQtq93qpysnNn+5WUPNWYFttceN5B3XFHeiVYKV7CflqP6rrKtMpakK69Vtq+vfPZRAAADFA3lZHJXh9MW/JWQjrdb1zlWLeubgpv25WUVvvMsc3tWpA+E8LJ3OQtWjr2C6pWSkzLmtb8vOtNepWu02dUvvbasD+CCQBgjA1HA2u/NZrC280wT3ofc3PVHpRGU4cTj7P4rsPy2TOi8+gkm15Ny1qree3S9L9/RevmX6GF/S+ubbgFAGCMrIwwkpaeStxuOMkKK+99b3X/je5z6FC4fvrpJ0PRol5z8tw5U5OPqxpKVunI46t1ULPaqhu08PDl4dfMxAEAjJnRCCOt1iFpV7tTiTNWaq0LK3HvTXKNkjaGhsplafm0J6mieZmO1xyOa0Lzx96pNWtUXykhnAAARp27F3K5QHKfmnJ3D382u14quU9MhD+zbp+ZCddnZtwrFXez6vWs7Wdnw/XZ2Xzbt7oe788s7K+dY8i4T2XqjT4z4z6lwy4d95BikpcTPq1HvFJJPRcAAAoi6U7vMBOMRmXk1lul48fDn2lxRWB2Nl9zanK4xL39IZk8y8q3cwzxkE3iPieXorfTVNFWTemR1JOuNruuPbhPC9pSrcxQKQEAjJpOU0y3l7YqI8nr6apDujLQqgKRvt7u9r2ogqSv57hPZdONPqMHfVJHXDpRXylZc9QrpT2hgiS5r18f9h3/CQBAH6mLysjohZFWYWLQ11uFkzxhpdV9KpVk8vBKaY9PTR7NDiU67BVtqQ09DOEAAPqsmzAynMM0yYbV9BBHekikn4/ViyGaPEM28cJoje6TnNZbKqm8eNHJZtcpHVb9tOCKLvzdn64djpIYwgEADKdOU0y3lwuaVQqSDau9rIRkVRzioZ48QzKdDNHkHbLpopJSKe3xKT2cqpSc8JJurh0KolICAOgTjewwTTw7Jjkbxr32Q7PbsJHVX5IOAo36TXo9RNPsejuBp8HzqGiLT+rRRCg54dJxn5720E+S3g89JQCAHhndMBKHgGZTd9sJCOvXh58nJvI1qzabItzqsdoICbmvZz1mO4En2qY0eUvrfpJSqXojza4AgC6NbhjJWvfDvXkYaRYYWq1H0s56Jq2qGung04v1SvI0v+YcOirpZm+0Rklp0321v6TZFQDQpW7CyHA0sMYNnY8+2vq8MMnVU9PbNluPJC2rkXT37nD/224LJ6dzD4+V1VR68cXhts2bO1tCPut55ml+Ta5PIjVcr2TRXiuv/LEqpb11Ta5LdzxVa3U4rE9SKtXux2l2BQAMWKcppttLw6m9edcBmZ1tbxgmb/9IfD1dZWnW15J+rOTtnTbDdlEVydpHRVsypgOHn0uTt9TvJ7leCQAALWjkh2k6HZbpJnzkbU7N6mvpplG1D0Gj5T4S9wnDN/X9JCXd3Hg/9JMAAFoY7zDSLAS008zazSycdKWg20bVAVdFsmbdhCqJe3LmzclVXBvtJw5mhBMAQMp4h5FmwzLNTpDXrJm13SGbbtY6Kboq0myKcEalxHQszLhptg4KQzgAgJTxCiPtDMs0Wxyt2Zl5+xku2g06RYWVRKioH7qJzgqcrJI0e2yqJACw4o1XGGlnzZG8i6P1csim11WUoqoi6RPxlfa46VhdL0nN2iTN9sOUYABY0boJI8MxtTd5fpi803OT01zd629LTptNT8Nttm2r6bXJc9FIzafaSq2n97a6f6vz1rQx3ffkfebmpIMHa+5T3v1C7dKVmee62apdWlCD6ctMCQYAdKvTFNPtpaYykqyGNKtSJIdwmvWP9GrIptdVkk6v5xkSaqeykix/ZNyn0blupOOtZ93QTwIAK5JGdpim1blomq390ax/pNMhm3bDRreNq3m37zSsJF/brPtkzVBK3KeiLW4Z04ArpT35+1LoJwGAFWF0w0irvpC8J9BLfyNv9iHfbKpwL8NGu30cPZy623Vja2KbSsV9Soc9LC2fmAYc95LEr2X8XmU9FlOCAWDs9T2MSHqxpPsk3S/pzRm3ny7pLyTdLeleSa9rtc+aE+U1CiN51xhJfqCWSvk+5Fut4Npt2Gi1eFs7j521fRdn923rPoltZnWgrkpS2nRfbbUlfr/icJh8/RjCAYCx1dcwImlC0pclfb+kNVHgOD+1za9L2hH9fKakA5LWNNvvBY0+eLN6P1qtMdLpMEyzdUraCRt5pvvG+8xTsehDhaOj+6S2Cb0kh71+GvDhxmuTzM7WhpPkiQWpkgDA2Oh3GPlRSR9PXH+LpLektnmLpHdJMklPiSooq5rt98k6w9drn5cmb/EZPeizejD0IjTq/Uh/+Cc/3PNWU5r1oAyikTWuBDXbPl3ZabV9v6oiTYZeKtri6TMCm46110tClQQAxkq/w8grJN2YuD4n6Y9S25wm6ZOSvi3psKSXtt7vBSe/WSe/ZYcPueTlhK/Ro16avMXX6xtemXpjvtkxeYZhkh+InVZJ0uEhT5Wkn1WRVv0b6SGTZKUi6ziyhl4mJkKVpO7Eew1OukejKwCMvX6HkZ/PCCP/O7XNKyT9r6gy8oOSvirpezL2NS/pznA5p+ZDrPXlROLPEFSmdTiEk0aBo52A0WmVpJfTeXtRFWkRIhruI7lN8jiaPe7ERDTj5ljN+9Ry+m+j58jCaQAwsroJI3kWPdsvaUPi+npJ30pt8zpJH4yO5/4ojDw9vSN3r7j7RnffOKF1mlx1TNIJrdZjmtRjkk5kXDy6tyX+XCVplZa1VvNHrpPpuMyPyea3yI4c1ho9qnVHvqmF834/LHZ26aWtFy/bvTsstnbbbWE7KXvhsmb7KZdr71sut7dIWvzYExOdLWh26aXV282k66+vLiIXX886hq1ba7e59dbqPnftavy4x4+rbH+sXaW/19Tk49F7ZVrSnNYuP6CFub31j3XVVY2fY/zasXAaAKwsrdKKpElJX1HoBYkbWJ+R2uZ6SW+Lfn6SpG9KOqPZfi+44IIQXZIl+oyfK6U9UWHiuK/S0ZrhmzzVlPg+J5ss81Q3klWSdENselZMuvrSrFG1VZUjvRR+t9OD2+klyTsM1GS/Wee4Kenm/BUZ+kkAYGRpAFN7L5H0JYVZNb8R/e5KSVdGP58l6a8lfV7SFySVWu3zZBjJKxVUKjO/FvUs1PeYNB/qOV7f19CslyTdEJsOG8nrzRpV8/SWJB83T69Jo6GVTmbY9Gj6b2nylrrXfHb6SO1J9/I0tjLrBgBGSt/DSD8ubYeRRjJCyszUEZ/Uo17bZ1IfSqTjYZ2MRtWFZL/F+vWNw0ZWI20/e0vyXM/T6NqqiTVPo2tGKKpMvbFuOfkZPVh7fKVS9Q1pFIyokgDAyFjZYaSRKJisn33YS1N/Fn041k5HrQsm5+5tHiA6bWTt9YybdmfY5B3S6WWjq+QzejC7KlWp1L4JjRZJywonVEkAYCgRRlpJfIiFL+SNhnSaBJNmH+atwka3M2xahYQ8wzHdDulkbdOs38TCImkzOuDpPpKKzdcGkWaLpGXtm1k3ADB0ugkjeWbTjL5yWdq3TyqXtbgoua+S+yqVSqb6GTthps7S154v82Oa081h1kejWSunn147+6TVtlJ71/PMhmk2Y+XSS5vPlsmatZOepbNzpzQ3Jx08mP24Wc971y6VL7pPB+17VdJi9DqH13jeb5DpuOY23SdddFF1PzMz0oED1ecsSbOz4TnEz0kKfy4shNvWrWPWDQCMuk5TTLeXgVZGWqiOwMSVkfRCbOH3J9c1GWRVpN11R/oxW6bZsEob/Saz00dSlSivXyem2fE1qpzEQzcM4QBAYURlpDvlsrS8LLlbRsXEFFdMwrom/0urDn4nVEwGURVpVsFoZ39Zx9don+nKylVXVV+sUim7WpNePyRahyS5zY4jV2tKh5WsRm098gda0JZ865A0qpxcfHG43/790rXX1r2/AIAh12mK6fYyTJWRRqoTPuLqSPJb/Qmfnng0rF1SRFWk370kPZrq22ib2irJ8TD9N6tK0umMH6okADBQooG1/5KfgZlnrS3taR0mmi1Pn2e2Szu35/nQ7uMCaK22qUy9MbWMvIeT7WlLe+uQcBI+ABgKhJEBib9wb9rkNR+i8WW1Hg1nH073lcQfklkLpGVVObL6UNqZDZP3fDTNFkDr1zokiW3CeW3S061PhLVfehmEqJIAQN8RRgqQXHak/nLCZ3Wg+i0//W29l8M5nVRN3DsLL33Ypjr993jN63dy+m83QSj5PJkODAB91U0YoYG1Q3HTa6UiTU2lbzUd1KzmVdHceXfUNnNu397bJlepvuG01TTddNNqq+nDiRPjdbWNVDvNePNmlT96mQ5qnSraqik9IkVNw/N+g9bqsBY2/239c2zQINvyJHycgA8AhlOnKabby6hXRtLiUZisSsm0HqmtkjRbPK3dptRhWQAtvU1yuff4cdILmSV7VKJtZmtWbXWfmTqSvRJrnuNp9FpSJQGAnhOVkeKVy2FNsEolzECdnKzetqxpzauitcsPaGHbvzauYtx2W/i2PzGRf4GzdqomWRWFPNOHm1UcTj9d2rtXmp+vfV5LS9UXIH6c9EJmyWnH0TY7dE00/dfD0zlySnXxuazHanbMjV5LD/umUgIAQ6LTFNPtZdwqI1mSxYFkP8Smyc9WGyuT39yTZwh2b69qkqeq0uzkdHln3GRVVpJPMD0rJz7bcZuVlUolvFY1fSSbbqx/rHQFJm8vC7NuAKCnRGVkOC0uZvWUmO44tlFr31TWwsXvldavr35z3707fJu/7bb2+zryLBufrFYkl1hv1ofRoudDu3dX91kq1d5nZiYs996qsiLVbVPWQrSMvJ983bbd8bL6x0o+p82bpW3bapetb9TLkqycbN9OlQQAitRpiun2shIqI0mVivvkZH2lZPXqaLbvzK/VVhv6NdU3rlZ0M+Om2aydHs+4qdi828kKSWJNl+R+0s8p+bokG3ny9JdQJQGAjoipvaMje+jG3XQiBJL08E0vVlptt6m12T7iBtB+Nb5mbFMp7fH0Crgnm4Kz9pO1iFujcNLoubM2CQC0pZswwjDNgGUP3Ugu09aHdmpB5TBsEA/fdNLU2myqb/r2rOGZrIbU5NDLjh1heCa5j04aX3NuU979wtSQTWgK3qYd2fvZtav2dUkOKcW/a3RG43hfnOcGAAan0xTT7WWlVkaSGk0Hnp5OfClPf0vvdVNrVpUkeTAFLoiW3qZS2uNTk0cTVZITPr3maL5z2rQ682/WKrRUSQAgNzFMM9qSn4ENQ0kseX6bfq7E2mlvyQC2qWhLzbBNrnPaNHqsVuGEXhIAyIUwMgYqlezl5ePP05oN42/q6RPvdbuYWS8CQ57+kx5sU5q8xdNTfzPPadPuY5nV98RQJQGAlggjY6RRKImXHqnbOP5wbDU80+pDutWMnbyBodnaJunHiY91Zqb6u6xtGgy91J/5N+OcNt0OF2UN4bCCKwDUIYyMobaGbtxbD9/0orekl/0n6WNNBpH0NlmBJnq+FW3xKR32ZB9JSTf3Zrio2RBOeoE6AFjhCCNjKvfQTbxx1vCN+2CGZ1qdiya9DkijoZCsbZJPvsGU3dpz2kRDNr0eLkoO4SQfHwBAGBl3jdYmmZ1tUCVJBpM8oaHdD+isfTQKDL3qLYkDTXrGSxR6QoXk4foKST9n/CSPKQ519JMAWKEIIytE1tBNZpUkKb2AWqNhk07PRZNVAemkt6TbbaJ0Vjp3b20gmbwl+/w1nYQj99bhJA6BNLoCWGEIIytIo16Shq0LyQ/GOJikh016PfW33d6S9ON0uU1JN9dXSJJVm25OrpdnCIfpwABWIMLICpNsC8k1bJO8Y6PhmzwfviM79TcKJHHVJp3i2j1/TbsBhioJgBWAMLJCZTW4thy2iSWHb7I+WHvRWzKI4Zkm26QDyfSao9Wpv+nKTr/CEVUSACtEN2GEc9OMsHJZWl6WZmerv3MPp5WZm2tx5+T5b6T6c8YsLVW3jc/r4t74fDZS63PepB+nXG69TbzfVufnydhm8fgViXPamJaPrtZWf5cWbD6cIOiii+rPX7NtWzhfTS+ea/KYtm8P593ZsCH8CQCo6jTFdHuhMtI7jYZtci+B0egswXl7S/o59NKDKklp031eu3z8Ca+U9jRvhk2tZ9LW83BvfkxUSQCMITFMA/f6dojcgSSryXVYhl561Oha0Zaa1VpNx8P5bLL2k7XmSRsrw+YOYvSSABgjhBGclNVH0nDV1kY7SDa5dhsYsqoLTVZUbWubNoNAxebddDwRSI6FM/622k/WTKJeNbpSJQEwJggjqJO1UFrba3K1WyXJCgxZ1YXkQTWqQKTLO8kg0E6VIhUEKqU9deezKU3ekn8/8QvZg3BElQTAOCGMIFNWIIm/0OfSTpUkT2BotUBaq7VMkkEgXaVoY2Gz9JDNyam/na4M20U4okoCYFwQRtBQ8nO07V6SpHSTa7MKSDow9KtvJPlY6WpLi4XNKqU9PjV51GvWItl0X0ervnYcjpLPJfl4VEkAjKBuwghTe8dcuSzt2iXNzNT+fmmpzRmmyanAt91WnRI7MSFdemn1+syMdOCAdOutvZuy616dQruwUPtYO3aEacexUqn6BGPxlN2DB0/uu3zRfVr2qZqpv0t3PFVzB6/Lf4zp5yqFedaXXlr7+I2mAyefS/Lx9u+Xrr02/Q4AwPjqNMV0e6EyMnhZwzZtV0jcq3OJG1VAip6B02yYJ2NIZXrN0erNerz9FVbbPblfO8+FKgmAESEqI8hjcTF8QiYXSVtayrFAWlq5HKoMBw6EKkm6ctDu4mfxfTpY2KzpNo0WNosXcIvud90rP60pHZZ0Qse1Smsnjmhh73n5Fmhr9lzjx4+3ife1d29YmS7Pc6FKAmAFIIysQDt2SFNT1etLS+Gzs+1QItWv5NrJ0Es3K5y2s82uXZlDKuXdL9SyTtOETkhapeWjqzW/9ALN+U2th1TyhiOz9oZwks+F1VsBjLtOSyrdXhimKV562GZiossdZk0F7sfCZt1MoW1y1t3ac9m4K16ptZ9DSG3MAmLGDYBhJoZp0InFxWq/p1T9It5RhUTqTZWkUQWi35WUXbu0+K7DqtiV0ZBNaGrduvR8Lcx/pn9DSFmNts3OjUOVBMA46jTFdHuhMjJckl/YpTZXbc3SqkqStQbJzEz1dwNYIr5ZJWV2+khNhaSkm9s7xl412qanUCefy+xsV+85APSSWGcE3ep6gbS0ZgumZa3LkeeDuN0hjS5m9zRdqbXZMeZ5fPfm4ajRuXEaPRdm3AAYAoQR9Ew6lHRdIXGvr5JkLRrWbBXWRkvE97JKkREWKtriU3rYaxZG082DO5FeO5UjqiQACkYYQU9VKrWf+11VSOIdJr+9x9WEVlWKZFjpZEijR42upU33ed1KrVlhoR8n0st6TZo9X6okAApCGEHP9aVCEkt+aHbb75H3vDhdBoGSbq6vkLQKUP0IR3meL1USAAUgjKAv0hWSuCjRU62qJO30ezSqpPSo36R26m8USFrtp0/hKNfzjR8XAAaAMIK+6dkS8o20qpL040R6XfSb1AWSTfd1FhZ6NIRUs41Z/fPtevEYAMiHMIK+qlTcp6b6GEhivaiStKpS9KDfJD1kU9GW9sNCr6skzaYMVyrVGUv0kwDoE8IIBqKvfSTuvamStFNJ6WJWTLqpdXrN0bBaa7thoZdTlt1bP178+tLoCqDHugkjrMCK3NIrti4vS1df3cMHKJelffvCn+7hd/04kZ579Vw1yW0anEhPUlgVNbGi6+Lr92hWB6MDNy0fXa1tSz/c/ISAyZPtSY3PVROfSC8+zrwryjY7f058bpyLL+YEfACGT6cpptsLlZHR1fcKifvgqiRdNMM2XIekmyGkdNdwJxWQdp4vVRIAPSIqIxikrArJNdf0+EEaVUniqkGvzxWT3mZmJpzeePfuhmfdLetGLVd2q1L6lBSdy2ZJc1q7/IAWrv5Ce+fGic9Vc9VV1degVMquknR7BuHkNlRJAAyDTlNMtxcqI6NvIBUS99pv7+lKQj/OFdNOc2j0ItQ2tbrP6MHsCki7PSHpruF+Pl+qJAC6IBpYUZT052X8Wdc3yfUzBrgKa56wUDp3b/2QTS+DT7+fb/y79ev7+AYCGFeEERSq72uRNBJ/EA9oFdY8YaGiLbWBJHlyvUGuwkqVBMCAdRNG6BlB1xYXQ8uDWfV3S0uhtaGvdu6U1q8Ps2JuvbV2lkrWrJh+9ZvE/R6HDqmsG1XSok72kBy7XAva0t6smFbb5H2+eZ5LcoZP+vHoJQEwKJ2mmG4vVEbGT6VSuzha34dskpILprm3XoW1Usle36NH/Rd1K7Weu7d/Q0iNnm8vzyBMlQRAC2KYBsMk+bk4sCGb5AdmJ1Nos/ovkk+m0fLzTT7gZ/Vg4iFSJ9fr9ZBKs+cb/67bx6OXBEAThBEMleTnXrIQMTDr1+f/gO22ObTJvis271OTR72uqbXRuWqywkI364skqyS9XPWVKgmADIQRDJ30kM3AKiTxg+epkrTTHNooeOQIC+lz2dQ0tRZ9Ir1kVaidKcNUSQCkEEYwtAoZsklqp0rSrLLQZb9H5sn1huVEet1MGaZKAiDSTRhhNg36ascOaWqqen1pSZqbG+ABbN8eZtx0uwpr/LuJibAKa+pcNSf35Z45K2ax8pim1xyLDsq0VTdo4dIP185k2bWr+nju1eMcwCygcFiJc/Pkfb7MuAHQC52mmG4vVEZWlqyRgYHrVZWkw/6Lis276djJTS2ukPR7SGVQz5cqCbCiiWEajILCA0k3vSQ9CgsVbXFLLBs/oce9YvODGVLp1ZTh5POtVGp/Ry8JsGIRRjAyCg8ksW6qJO2EhYwP+UppT83Zfk0nvFLa099G235MGW42RXp2lioJsMIQRjBSBnaCvWZ6WSVpFhYaVEnW6xs1r4HpWBiy6XWjrXvHs4C6DmPx68vwDbAiEEYwctKBZKCrtab1ukqSFRZSQyoVbfEZHfC6M/22ExaaDan0ewin3coNwzfA2COMYCQNTSDptkri3jwsxL/LCAuhhyRuaj3h02uOhh6SXgaBZDgZxJThZttQJQHGFmEEIyv5GZr8bC9MJ1WSLsNC7Zl+3dfrG8XOiul1lSS9DVUSYCwRRjDS0oGkkB6S5ME0qpK4920KbencvS4dP3kpTd7SnyBQ9JRhqiTA2CKMYOSlJ2UU2kMSS1ZJBtB/UdvUmjix3rCvwtpNf8vsbHHvL4Ce6iaMsAIrhkK5LJVK1evu0rZtxR2PpNrVW2+7retVWFtts13XakqPSHJJpiXNaU43D34V1uuvl269tbPnu7CQ/XiNtnEP1zdsCH8CWJk6TTHdXqiMIEt6yKawdUjSKpXq+hl97r+oO7HeuXtry0b9HlJp1Qzb7CzD3VRuqJIAI01URjAuyuVwipbY0pK0du0QfGkul6WDB6UDB0KVJK4QzM6GqsHSUnXbLqski3qNSlrUyQrJ154fKiRSKB/FVYv4fp1WScrl1se0a1d4vOPHw23x8423iffVi8oNVRJg5eo0xXR7oTKCZtL9lYXPsklKN7mm+y96VCUpbbqvtkKy6b7sikTW4/VjFtAgm2GpkgAjRzSwYhwNzdLxzaSnAqc7cbv8YE4P2UzrcO1KrY0eL53i+jkrpp/NsMy6AUZGN2GEYRoMrcXFMCoQW1qS1q0bsgp+sslVCsMQsawhlc2bQ2fuwYNhm2bDJYcOhSGbyfdFOzQta622aUfrx5Oqjaezs7WPF99vFJphr7pK2r9fuvbanG8IgJGUJ7FIerGk+yTdL+nNDba5WNLnJN0raU+rfVIZQR6VivvUVO2X/6EasknKqpK0WoU153BJ5pDNKKzC2kkzbKNqElUSYKipn8M0kiYkfVnS90taI+luSeentpmR9M+SzomuP7HVfgkjaEfyM3xoh2wa9ZI0+2BuYwindqXWaB2SYVmF1b13K8M2GuaJH48VXIGh1O8w8qOSPp64/hZJb0lt8wZJv9POAxNG0I54Zu3QB5JYVpWk2Qdzzt6KugpJJyu1Nmo8HZZm2DzVnTi0UCUBhka/w8grJN2YuD4n6Y9S21wn6Z2SbpN0l6RXt9ovYQSdyPr8HEqdVElyDqnUrUMyeUv+ikSzxtNhaobNW0mJX2OGb4DCdRNGJnO0lVhWq0nq+qSkCyT9lKQpSZ82s9vd/Us1OzKblzQvSeecc06OhwZqLS6GP+NlPebnpTe9SbruutAHOjTK5eoBbdhQ3wiabNa8/vpqA6sUGk5f+lJp9+6wjVRdy2NpSdFLoCXNSTItHbtc0uNaPP3q6r5brS+ybl21gTRe2CVvM2y3ja6dbDMzE9Z4iY87fp0uvrh6v2uvHbK/BADyyjObZr+kDYnr6yV9K2Obv3L3R9z9u5L2Snp2ekfuXnH3je6+8cwzz+z0mLHCpWfZLC8P+WSLrBk3zT50d+wIQaTJwmqLlcdU2vSvUnLp+PPuyLewWqMl25vNitm8ORxTn5fEz33c6ddpYiK8zgBGU6vSiULV4yuSnqJqA+szUtv8kKS/ibadlvQFSc9stl+GadCN9CybQs/02452eklaDamY1Q3ZVLSlv0Mq3TbDNjqmXpwdGUCh1M91Rtz9mKRfkvRxSV+U9Kfufq+ZXWlmV0bbfFHSX0m6R9JnFHpMvtDDzATUKJdDRWT9+nB9eTl8SR+qNUiytFMlyVqDJK5aRFWD6jokoUIyr0pYOj5PBaTT4ZIclRtdf311ifg8VZl2KjeNXicAo6vTFNPthcoIeiH5pXrkviTnqZLknIEyqwcTvadRU2u/TqTXTTNsv44JQOHEcvBYydKBZOin/cbSs0C6+GCuaItP6eGTQzamY9X7dbocey/WF4nfjDzTmDtd8wTAUCCMYMUb2UCSFH84dxEWptccrVZHdHP/V2FtN8A0m8bcSSUFwNAgjABe39Q6coEkWSnpsEpSsXnPXIOkWeNpOiz0sqm0k2bYvJUUAEOlmzDCifIwNuKm1lKp+rulJWlurrhjaku5LO3bF/50D79rs/G07BWVtKi4oXXp2OWhoTVep+PWW2vX6bj00moDabyvdhpPe7F2SHxMyQbdVifb27mz01cZwDDqNMV0e6Eygn7KqvqPlC6rJJkn1svZDDvQxtNOhnkYngGGkhimAeqNfCCJpXtJcg6p1C0br5vzh4V2hkt60ejKzBlg5HUTRhimwdhaXBzhIZuknTvD2iQ7d4alZj1jCMfrh1TSa5AsaS4M2eQZUun1Kqy9GuZhPRFgPHWaYrq9UBnBoIxNhcQ9VAbiM9bmnBVTmrylpkJSKe3p3yqsjSo3vVjPBMBQE5URoLGsCsm6dSOwWmuWclk6eDA0ft52W3bjaarasHj8ipqm1vmlF2ju+Lvzr/raqvFUat0M2+2J9KiIAOOt0xTT7YXKCAYtXSEZ+dmh6SpJi8bT2ekj1fYPPd67ptJummHzrnkCYOiJygjQWvpsv4cOSWvXjmiFRKqvksRKpepUWI/6NjZv1g6/RlM6LOmEjmuVbL6suTec2vuz7u7aVT1/TYtj0u7d9T0w6R4UqiLA2COMYEXZsSOMMsSWl6Vt2wo7nN6JT8DXZJ2O8pF3aFmnaUInFP/TXzp2uRa0pXdrh+Rtho3XDmk1zEMQAVYEwghWlLiYkOwheeihEa6OxJILpsUzbiYmQpVCqgkCm8/9tKQTintItuoGLZz3+/lmxbRbJcnaJn1Mjc4EzMJmwIpBGMGKtLgYvrCbhc/IrVvHIJDE4ipJgyCw+NBlck1oVgclSa4Jbb3jtVrw14fte9142skwz4EDVEWAFcQ8LokO2MaNG/3OO+8s5LGB2Lp1oVIihS/o8c9jY8MGaf/+2rCwdavkroWpN2rrkevkMkmS6bh2Tf2Kyj93IEw5ilUqIZzEv4v3dc014QWLA4SUHU6ix9PMTAgerbYZuzcBWBnM7C5339jJfamMYEXbsaP686FDI7ooWjNNqiTlP3ymdk29SabQt+Ga0DZlNJ5KteFk8+bQaBOHhk4rII22AbDiEEawopXLtTNslpbGaLhGqu8lSVUkyo/+b+3SlTKdkCQdOnKK1h5/KDS1xo2nV11V3V9yVoxUXYV1drY2nEissAogN8IIVrwdO6Spqer1+fkxrJBIDask5Yk/0a7S3yluaF3WWm3VLi3sPa83S8THlZR4KKZRlYQgAqxY9IwAkWT/iBQ+e8f28zGjl2Ru/hQtaU462UNyQru0VWX748Y9IXHfiBSqIy99aXVWjBR6QMyq28Tri3z0o7WzaegTAUYePSNAD2RVSEZ6UbRmMqoki3qNKnalTOELimtVmPb7hF8+uU3DYZdmS8RnrS+SrKTQJwKseFRGgJR0hWT9+tB2MbZSVZKFvedp69KPyzUhSZqZelQHn3BW72fOxPsa2/ITsLJQGQF6KF0hOXBgTKsjsVSVpLz7hVFTa6hwHDpyiuYOXhe27eXMGYIIgAhhBEgpl8My8evXh+vLy2Pc1CplzrgJTa1/r7ipdUlzWqvDWrj0wyGZdbtEPA2rABIII0AD27eHL/WxsZv2myVRJSkvXqTS5PuUnGWz7QPP680S8QQRAAmEEaCBcjmMJCSHbMZq2fgsySqJpMVT36CSFqWoqfXQkVM0d/zd3S8RDwAJhBGgiXjIJl4YzX3Mh2zSdu7U4vrfUKX0KSWHbOY2fLJ+7ZD0+iKssAogJ8IIkMOOHStwyEY6WSkpL16k0rnVhdGWvvZ8rT20P5xcL147ZPfuEDqkapUkHU6oigDIQBgBcliRQzYpi199gUqlOJFFPSTaUV07JL2+SFwRkRieAdAUYQTIKWvIZsUFksX43HkhZDy2+rRwQ3Lxsx07qquwmoUXjOEZAE0QRoA2JYdsVmogmZ0NL8CRx1dr7fwV1eGa9MyZXbvCQi1URQA0QRgB2hQP2azkQLJjR/XnZU2HE+tt/tvwC2bOAGgTYQToQFYg2bat2GMapHI5Hq4JXKu07aMXhUXTmDkDoE2EEaBD6UBy6NAKmvKrMFxTqaSe/9l/U11anqoIgJwII0AX4kASW1paWYGk7vnf8TTZ/n2a20sQAZAfYQToUnrIYiUGkk2ban+3e3cxxwJgNBFGgB6oTnkNlpaktWtXRlPr3Jx0xx21v9u8uZhjATCaCCNAj6QDyfJy6OccV3NzoV9kaan6O7PQR7K4WNxxARg9hBGgh9KBZP/+8RyyWVioDSGSND0d+kfoWwXQLsII0GOLi2FCSWzcekjm5sLJApNKJemRRwgiADpDGAH6YPv22vPYjPKJ9RYWpA0bQgiZnq4flnFnWAZAdwgjQB/E57FJDtmM6iqt114bhpuWlqQjR2pvu+KKYo4JwHghjAB9FM7jEn52D8MbozZkc/HF9b8rlaiIAOgdwgjQZ8kT60nSe99b3LG0Y2GhflhmeprZMgB6jzAC9Fl6lVL34V6DJA4h8/O1wzKVCk2qAPpjsugDAFaC+AM8noWyvBx6SJK3DYOFhXBc7rW/L5WG6zgBjBcqI8CA1J3pdsjO9JsVRBiWATAIhBFggNJnun3ooeEYronXDomDSLySKsMyAAaBMAIMWNxDEq/RUeQMm6wmVTNWUgUwWIQRoADlsjQzU71exCqtcTUk2aTKku4AikAYAQqyY0f9Kq2DCCRZ1RCJJd0BFIcwAhQka5XWfgeSRtUQmlQBFImpvUDB4hAQVyriP3sZDhYWpKuvrl/OvVQihAAoHpURYAgsLvavQhJP2aUaAmBYEUaAIdHLQBKfaffCC2un7Er0hgAYPoQRYIhkBZJO1iGJz7R7xx3V38Vrh1ANATBsCCPAkEkHkq1b2wskc3MhiCQxZRfAMCOMAENocVGanQ0/t7Mw2oUXMmUXwOghjABDaseO6rLxUvMeknjtkOSwDE2qAEYFU3uBIRVXMpJTcrOm/Wad4G7TJun22wdznADQLSojwBBrtDBa3EOSPsGdFLYliAAYJYQRYASkm1rn58MQTvoEdwzLABhFhBFgRCSbWtOYLQNglBFGgBGSPrmexGwZAKOPBlZghJTLhA4A44fKCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChcoURM3uxmd1nZveb2ZubbPdcMztuZq/o3SECAIBx1jKMmNmEpHdKeomk8yVtNrPzG2y3Q9LHe32QAABgfOWpjDxP0v3u/hV3PyrpfZIuy9julyV9QNIDPTw+AAAw5vKEkbMl7Utc3x/97iQzO1vSyyTd0LtDAwAAK0GeMGIZv/PU9eskXePux5vuyGzezO40szu/853v5DxEAAAwziZzbLNf0obE9fWSvpXaZqOk95mZJJ0h6RIzO+buH0pu5O4VSRVJ2rhxYzrQAACAFShPGPmspKea2VMkfVPS5ZJeldzA3Z8S/2xmN0n6aDqIAAAAZGkZRtz9mJn9ksIsmQlJ73b3e83syuh2+kQAAEDH8lRG5O4fk/Sx1O8yQ4i7v7b7wwIAACsFK7ACAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABQqVxgxsxeb2X1mdr+ZvTnj9ivM7J7o8g9m9uzeHyoAABhHLcOImU1Ieqekl0g6X9JmMzs/tdlXJf2Euz9L0tslVXp9oAAAYDzlqYw8T9L97v4Vdz8q6X2SLktu4O7/4O4Ho6u3S1rf28MEAADjKk8YOVvSvsT1/dHvGnm9pL/s5qAAAMDKMZljG8v4nWduaPaTCmHk+Q1un5c0L0nnnHNOzkMEAADjLE9lZL+kDYnr6yV9K72RmT1L0o2SLnP3B7N25O4Vd9/o7hvPPPPMTo4XAACMmTxh5LOSnmpmTzGzNZIul/SR5AZmdo6kD0qac/cv9f4wAQDAuGo5TOPux8zslyR9XNKEpHe7+71mdmV0+w2Stkv6XknvMjNJOubuG/t32AAAYFyYe2b7R99t3LjR77zzzkIeGwAA9JaZ3dVpIYIVWAEAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACkUYAQAAhSKMAACAQhFGAABAoQgjAACgUIQRAABQKMIIAAAoFGEEAAAUijACAAAKRRgBAACFIowAAIBCEUYAAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEIRRgAAQKEIIwAAoFCEEQAAUCjCCAAAKBRhBAAAFIowAgAACpUrjJjZi83sPjO738zenHG7mdk7otvvMbPn9P5QAQDAOGoZRsxsQtI7Jb1E0vmSNpvZ+anNXiLpqdFlXtL1PT5OAAAwpvJURp4n6X53/4q7H5X0PkmXpba5TNJ7PLhd0oyZfV+PjxUAAIyhPGHkbEn7Etf3R79rdxsAAIA6kzm2sYzfeQfbyMzmFYZxJOkxM/tCjsfH4Jwh6btFHwRO4v0YPrwnw4X3Y7ic1+kd84SR/ZI2JK6vl/StDraRu1ckVSTJzO50941tHS36ivdkuPB+DB/ek+HC+zFczOzOTu+bZ5jms5KeamZPMbM1ki6X9JHUNh+R9OpoVs2Fkh5y9293elAAAGDlaFkZcfdjZvZLkj4uaULSu939XjO7Mrr9Bkkfk3SJpPslLUt6Xf8OGQAAjJM8wzRy948pBI7k725I/OySfrHNx660uT36j/dkuPB+DB/ek+HC+zFcOn4/LOQIAACAYrAcPAAAKFTfwwhLyQ+XHO/HFdH7cI+Z/YOZPbuI41xJWr0nie2ea2bHzewVgzy+lSbP+2FmF5vZ58zsXjPbM+hjXGly/L91upn9hZndHb0n9C32kZm928weaLQ8R0ef6+7et4tCw+uXJX2/pDWS7pZ0fmqbSyT9pcJaJRdKuqOfx7SSLznfjx+TNBv9/BLej+Lfk8R2f6vQu/WKoo97XC85/43MSPpnSedE159Y9HGP8yXne/LrknZEP58p6YCkNUUf+7heJF0k6TmSvtDg9rY/1/tdGWEp+eHS8v1w939w94PR1dsV1oxB/+T5NyJJvyzpA5IeGOTBrUB53o9XSfqgu39Dktyd96S/8rwnLuk0MzNJpyqEkWODPcyVw933KrzGjbT9ud7vMMJS8sOl3df69QrpFv3T8j0xs7MlvUzSDUK/5fk38jRJs2Z2m5ndZWavHtjRrUx53pM/kvRDCottfl7S1e5+YjCHhwxtf67nmtrbhZ4tJY+eyP1am9lPKoSR5/f1iJDnPblO0jXufjx88UMf5Xk/JiVdIOmnJE1J+rSZ3e7uX+r3wa1Qed6TF0n6nKQXSvoBSZ8ws0+5+3/0+diQre3P9X6HkZ4tJY+eyPVam9mzJN0o6SXu/uCAjm2lyvOebJT0viiInCHpEjM75u4fGsgRrix5/8/6rrs/IukRM9sr6dmSCCP9kec9eZ2k3/PQsHC/mX1V0tMlfWYwh4iUtj/X+z1Mw1Lyw6Xl+2Fm50j6oKQ5vukNRMv3xN2f4u7nuvu5kv5c0hsIIn2T5/+sD0t6gZlNmtm0pE2Svjjg41xJ8rwn31CoVMnMnqRwwravDPQokdT253pfKyPOUvJDJef7sV3S90p6V/RN/JhzIqq+yfmeYEDyvB/u/kUz+ytJ90g6IelGd+cM5H2S89/I2yXdZGafVxgiuMbdOZtvn5jZbkkXSzrDzPZL+i1Jq6XOP9dZgRUAABSKFVgBAEChCCMAAKBQhBEAAFAowggAACgUYQQAABSKMAIAAApFGAEAAIUijAAAgEL9P76n1BPpD3CsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(9, 8))\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_ylim(0,1)\n",
    "ax.set_title(\"MPM (blue) vs FEM (red) beam bending, E=500, nu=0.35\")\n",
    "\n",
    "ax.scatter(points_FEM_X, points_FEM_Y, color='red', s=3)\n",
    "ax.scatter(points_MPM_X, points_MPM_Y, color='blue', s=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab66d5a3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bd10636",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa5c57ba",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

MPM_beam_bending.py

import numpy as np

import taichi as ti

ti.init(arch=ti.cpu)  # Try to run on GPU

quality = 1  # Use a larger value for higher-res simulations
n_particles, n_grid = 100000 * quality**2, 128 * quality
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 0.1e-4 / quality   #1e-4
elapsed_time =ti.field(dtype=float, shape=())
n_substeps =ti.field(dtype=int, shape=())
indice_temoin =ti.field(dtype=int, shape=())
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho
E, nu = 500, 0.35  # Young's modulus and Poisson's ratio  0.1e5, 0.2

mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / (
    (1 + nu) * (1 - 2 * nu))  # Lame parameters
x = ti.Vector.field(2, dtype=float, shape=n_particles)  # position
v = ti.Vector.field(2, dtype=float, shape=n_particles)  # velocity
C = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # affine velocity field
F = ti.Matrix.field(2, 2, dtype=float,
                    shape=n_particles)  # deformation gradient                 
material = ti.field(dtype=int, shape=n_particles)  # material id
Jp = ti.field(dtype=float, shape=n_particles)  # plastic deformation
grid_v = ti.Vector.field(2, dtype=float,
                         shape=(n_grid, n_grid))  # grid node momentum/velocity

grid_m = ti.field(dtype=float, shape=(n_grid, n_grid))  # grid node mass

# geometrical dimensions of beam 
length_beam = 0.6; 
width_beam = 0.10
point_bas = 0.8  
indice_bas = int(point_bas * n_grid)
indice_haut = int(indice_bas + width_beam * n_grid)
gamma = 0.9999  # damping coefficient
g = 15 

@ti.kernel
def substep():
    for i, j in grid_m:
        grid_v[i, j] = [0, 0]
        grid_m[i, j] = 0
    for p in x:  # Particle state update and scatter to grid (P2G)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        # Quadratic kernels  [http://mpm.graphics   Eqn. 123, with x=fx, fx-1,fx-2]
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
        F[p] = (ti.Matrix.identity(float, 2) +
                dt * C[p]) @ F[p]  # deformation gradient update            
        h = ti.exp( 10 * (1.0 - Jp[p]))  # Hardening coefficient: snow gets harder when compressed
        if material[p] == 1:  # jelly, make it softer
            h = 0.3
        h = 1.0 # no change
        mu, la = mu_0 * h, lambda_0 * h       
        if material[p] == 0:  # liquid
            mu = 0.0
        U, sig, V = ti.svd(F[p])
        J = 1.0
        for d in ti.static(range(2)):
            new_sig = sig[d, d]
            if material[p] == 2:  # Snow
                new_sig = min(max(sig[d, d], 1 - 2.5e-2),
                              1 + 4.5e-3)  # Plasticity
            Jp[p] *= sig[d, d] / new_sig
            sig[d, d] = new_sig
            J *= new_sig
        if material[p] == 0:  # Reset deformation gradient to avoid numerical instability
            F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J)
        elif material[p] == 2:
            F[p] = U @ sig @ V.transpose(
            )  # Reconstruct elastic deformation gradient after plasticity
        stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose(
        ) + ti.Matrix.identity(float, 2) * la * J * (J - 1)                     
        stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
        affine = stress + p_mass * C[p]                                         
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # Loop over 3x3 grid node neighborhood
            offset = ti.Vector([i, j])
            dpos = (offset.cast(float) - fx) * dx
            weight = w[i][0] * w[j][1]
            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)   
            grid_m[base + offset] += weight * p_mass

    for i, j in grid_m:        
        if grid_m[i, j] > 0:  # No need for epsilon here
            grid_v[i,
                   j] = (1 / grid_m[i, j]) * grid_v[i,
                                                    j]  # Momentum to velocity
            grid_v[i, j][1] -= dt * g  # gravity 50
            if i < 3 and grid_v[i, j][0] < 0:
                grid_v[i, j][0] = 0  # Boundary conditions
            if i > n_grid - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0
            if j < 3 and grid_v[i, j][1] < 0: grid_v[i, j][1] = 0
            if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0
            
            # definition of beam
            if i < 3 and (j >= indice_bas and j <= indice_haut):
                grid_v[i, j][0] = 0
                grid_v[i, j][1] = 0
    for p in x:  # grid to particle (G2P)
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[p] * inv_dx - base.cast(float)
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1.0)**2, 0.5 * (fx - 0.5)**2]
        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # loop over 3x3 grid node neighborhood
            dpos = ti.Vector([i, j]).cast(float) - fx
            g_v = grid_v[base + ti.Vector([i, j])]
            weight = w[i][0] * w[j][1]
            new_v += weight * g_v * gamma
            new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
        v[p], C[p] = new_v, new_C
        x[p] += dt * v[p]  # advection
    old_elapsed_time = elapsed_time[None]
    elapsed_time[None] += dt
    n_substeps[None] += 1
    if int(elapsed_time[None]/0.1) > int(old_elapsed_time/0.1) :
        print('elapsed time = ',elapsed_time[None])
    

@ti.kernel
def initialize():
    for i in range(n_particles):
        x[i]=[ti.random()*length_beam, ti.random()*width_beam + point_bas]
        if x[i][0] > 0.98 * length_beam and indice_temoin[None] == 0: # for printing evolution of one particle close to the end of beam
            indice_temoin[None] = i
        material[i] = 1  # 0: fluid 1: jelly 2: snow        
        v[i] = ti.Matrix([0, 0])
        F[i] = ti.Matrix([[1, 0], [0, 1]])
        Jp[i] = 1
    elapsed_time[None] = 0.0
    n_substeps[None] = 0


initialize()
print("info: oscillations will be small enough after elapsed time > 0.5; you can stop before by typing escape on beam image")
gui = ti.GUI("Taichi MLS-MPM-99", res=1024, background_color=0x112F41)
limit_range = 30 #int(max(30,2e-3 // dt))

while (not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT)) and  (elapsed_time[None] < 0.6):    
    for s in range(limit_range):   #range(int(2e-3 // dt))
        substep()
        if n_substeps[None]%400 == 0 :
            print("nb_steps = ", n_substeps[None], ", particule_temoin X = ", x[indice_temoin[None]][0], ", particule_temoin Y = ",x[indice_temoin[None]][1] )
    gui.circles(x.to_numpy(), radius=1.5, color=0x068587)
    gui.show(
    )  # Change to gui.show(f'{frame:06d}.png') to write images to disk

# opening file to write contour of beam
fich = open("fichier_plot_beam_MPM.txt", "w" )
scatter_plot_MPM = []
liste_p = []

n_pas = 1024
pas = 1/n_pas
for p in range(n_particles) :
    indice_i = int(x[p][0]/pas)
    liste_p.append([indice_i, x[p][1]])

# select contour particles
for i in range(n_pas):
    ligne_v = [el[1] for el in liste_p if el[0]==i ]
    if len(ligne_v) > 0 :
        min_y = min(ligne_v )
        max_y = max(ligne_v )
        scatter_plot_MPM.extend([[i*pas, min_y], [i*pas, max_y]]) 

# write contour of beam to file
for p in scatter_plot_MPM:
    fich.write(str(p[0]) + " " + str(p[1]) + "\n")
  
fich.close()