sloshing_quadrilateral.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e09cfced-6a1d-425c-bb38-bab6f6417630",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "def he_converge():\n",
    "    he = []\n",
    "    e1 = []\n",
    "    e2 = []\n",
    "    e3 = []\n",
    "    for i in range(6):\n",
    "        h = 1./2**i\n",
    "        he.append(h)\n",
    "        [a,b,c] = main(h)\n",
    "        e1.append(a)\n",
    "        e2.append(b)\n",
    "        e3.append(c)\n",
    "        \n",
    "    # converge speed\n",
    "    slope1 = (np.log(e1[-1])-np.log(e1[0]))/(np.log(he[-1])-np.log(he[0]))\n",
    "    slope2 = (np.log(e2[-1])-np.log(e2[0]))/(np.log(he[-1])-np.log(he[0]))\n",
    "    \n",
    "    # plot dx vs error in loglog\n",
    "    plt.loglog(he,e1,'-o',he,e2,'-x',he,e3,'-s')\n",
    "    plt.xlabel(\"he\")\n",
    "    plt.ylabel(\"RMS error\")\n",
    "    plt.legend(['potential','u','v'])\n",
    "    plt.text(0.1,0.1,'slope phi = '+str(slope1))\n",
    "    plt.text(0.1,0.05,'slope vel = '+str(slope2))\n",
    "    \n",
    "def check_ita(ita,h):\n",
    "    if np.abs(ita)>10.*h:\n",
    "        sys.exit(\"wave too high\")\n",
    "    return\n",
    "\n",
    "def exa_sol(af,L,n,x):\n",
    "    exa_phi = np.zeros(n)\n",
    "    exa_u = np.zeros((n,2))\n",
    "    for i in range(n):\n",
    "        exa_phi[i] = af*np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*(x[i][1])/L)\n",
    "        exa_u[i][0] = af*(math.pi/L)*np.cos(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*(x[i][1])/L)\n",
    "        exa_u[i][1] = af*(math.pi/L)*np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.sinh(math.pi*(x[i][1])/L)\n",
    "    return exa_phi,exa_u\n",
    "    \n",
    "def rmse(n,exa,num):\n",
    "    e = 0.\n",
    "    for i in range(n):\n",
    "        e += (exa[i]-num[i])**2\n",
    "    e = math.sqrt(e/n)\n",
    "    return e\n",
    "\n",
    "def contour_plot(nx,ny,x,sol,exa):\n",
    "    X = np.zeros(nx)\n",
    "    Y = np.zeros(ny)\n",
    "    Z_sol = np.zeros((ny,nx))\n",
    "    Z_exa = np.zeros((ny,nx))\n",
    "    for i in range(nx):\n",
    "        X[i] = x[i][0]\n",
    "    for i in range(ny):\n",
    "        Y[i] = x[nx*i][1]\n",
    "    c = 0\n",
    "    for i in range(ny):\n",
    "        for j in range(nx):\n",
    "            Z_sol[i][j] = sol[c]\n",
    "            Z_exa[i][j] = exa[c]\n",
    "            c += 1\n",
    "            \n",
    "    plt.contour(X, Y, Z_sol, 20, colors='black')\n",
    "    plt.contour(X, Y, Z_exa, 20, colors='red')\n",
    "    #plt.legend(['numerical reults','exact function'])\n",
    "    \n",
    "def plotMesh(elements, nodes,l2g=False):\n",
    "    from matplotlib import pyplot as plt\n",
    "    plt.figure(figsize=(18,9)); plt.rcParams.update({'font.size': 22})\n",
    "    for eN, e in enumerate(elements):\n",
    "        element_nodes = nodes[e[0]], nodes[e[1]], nodes[e[2]], nodes[e[3]]\n",
    "        X = [n[0] for n in element_nodes]; Y = [n[1] for n in element_nodes]\n",
    "        cx = sum(X)/4.; cy = sum(Y)/4.\n",
    "        X.append(element_nodes[0][0]); Y.append(element_nodes[0][1])\n",
    "        plt.plot(X, Y); plt.text(cx, cy, \"$T_{{{0:d}}}$\".format(eN))\n",
    "        for n in e:\n",
    "            plt.text(nodes[n,0], nodes[n,1], \"$N_{{{0:d}}}$\".format(n))\n",
    "    #if l2g:\n",
    "        #plt.text(nodes[n,0], nodes[n,1], \"$N_{{ {0:d} = elements[{1:d},{2:d}] = elements[0,1]}}$\".format(n,eN,2))\n",
    "    plt.axis('equal')\n",
    "    return plt\n",
    "\n",
    "def mesh(L,h,he,method):\n",
    "    # element number and size\n",
    "    nex = int(math.ceil(L/he))\n",
    "    ney = int(math.ceil(h/he))\n",
    "    ne = nex*ney\n",
    "    dx = L/nex\n",
    "    dy = h/ney\n",
    "    \n",
    "    # node number\n",
    "    nx = nex+1\n",
    "    ny = ney+1\n",
    "    n = nx*ny\n",
    "    \n",
    "    # mesh method\n",
    "    if method == 'even':\n",
    "        xi = np.linspace(0.,L,nx)\n",
    "        yi = np.linspace(0.,h,ny)\n",
    "    elif method =='cos':\n",
    "        xi = np.linspace(0.,L,nx)\n",
    "        yi = x_cos(0.,h,ny)  \n",
    "    elif method =='pwr':\n",
    "        xi = np.linspace(0.,L,nx)\n",
    "        yi = x_power(0.,h,ny)  \n",
    "        \n",
    "    # mesh node\n",
    "    c = 0\n",
    "    x = np.zeros((n,2))\n",
    "    for j in range(ny):\n",
    "        for i in range(nx):\n",
    "            x[c][0] = xi[i]\n",
    "            x[c][1] = yi[j]\n",
    "            c += 1\n",
    "            \n",
    "    # local-global index for triangle mesh (left-down right-up)\n",
    "    c = 0\n",
    "    ln = np.zeros((ne,4), dtype=int)\n",
    "    for j in range(ney):\n",
    "        for i in range(nex):\n",
    "            ln[c][0] = nx*j+i\n",
    "            ln[c][1] = nx*j+i+1\n",
    "            ln[c][2] = nx*(j+1)+i+1\n",
    "            ln[c][3] = nx*(j+1)+i\n",
    "            c += 1\n",
    "   \n",
    "    # define global index of the boundary node for solving BVP\n",
    "    epi = 1.e-6\n",
    "    idir = []\n",
    "    ineu1 = []\n",
    "    ineu2 = []\n",
    "    ineu3 = []\n",
    "    for i in range(n):\n",
    "        if abs(x[i][1]-h) < epi:   # free surface\n",
    "            idir.append(i)           \n",
    "        elif abs(x[i][0] - 0) < epi: # left wall\n",
    "            ineu1.append(i)\n",
    "        elif abs(x[i][0] - L) < epi: # right wall\n",
    "            ineu3.append(i)\n",
    "        elif abs(x[i][1] - 0) < epi: # bottom\n",
    "            ineu2.append(i)\n",
    "              \n",
    "    ndir = len(idir)\n",
    "    nneu1 = len(ineu1)\n",
    "    nneu2 = len(ineu2)\n",
    "    nneu3 = len(ineu3)\n",
    "    \n",
    "    # define global index of the free-surface nodes for time marching\n",
    "    ifs = []\n",
    "    for i in range(n):\n",
    "        if abs(x[i][1]-h) <epi:   # free surface\n",
    "            ifs.append(i)  \n",
    "    nfs = len(ifs)\n",
    "    \n",
    "    return nex,ney,ne,nx,ny,n,ln,ndir,idir,nneu1,ineu1,nneu2,ineu2,nneu3,ineu3,ifs,nfs,x\n",
    "\n",
    "def x_cos(x0,xn,n):\n",
    "    x = np.zeros(n)\n",
    "    r = xn-x0\n",
    "    ds = 0.5*np.pi/(n-1)\n",
    "    for i in range(n):\n",
    "        x[i] = x0+r*np.sin(ds*i)\n",
    "    return x\n",
    "\n",
    "def x_power(x0,xn,n):\n",
    "    w = np.power(np.linspace(0,1,n),0.7)\n",
    "    s = np.sum(w)\n",
    "    l = xn-x0\n",
    "    x = np.zeros(n)\n",
    "    x[0] = x0\n",
    "    for i in range(n-1):\n",
    "        x[i+1] = x[i]+l*w[n-1-i]/s\n",
    "    return x\n",
    "\n",
    "def load_phi(n,ndir,idir,phi):\n",
    "    F = np.zeros(n)\n",
    "    for i in idir:\n",
    "        F[i] = phi[i]\n",
    "    return F\n",
    "\n",
    "def loc_nod():  \n",
    "    # local coordinate of corner nodes\n",
    "    xc = np.zeros((4,2))\n",
    "    xc[0,0] = -1.\n",
    "    xc[1,0] = 1. \n",
    "    xc[2,0] = 1. \n",
    "    xc[3,0] = -1. \n",
    "    xc[0,1] = -1.\n",
    "    xc[1,1] = -1. \n",
    "    xc[2,1] = 1. \n",
    "    xc[3,1] = 1. \n",
    "    return xc\n",
    "\n",
    "def Shapefunctions(xi,eta):\n",
    "    sha = np.zeros(4)\n",
    "    sha[0] = 0.25*(1-xi)*(1-eta)\n",
    "    sha[1] = 0.25*(1+xi)*(1-eta)\n",
    "    sha[2] = 0.25*(1+xi)*(1+eta)\n",
    "    sha[3] = 0.25*(1-xi)*(1+eta)\n",
    "    return sha\n",
    "\n",
    "def dNdrs(xi,eta):\n",
    "    # ND = NaturalDerivatives = dSh/dxi and dSh/dita\n",
    "    ND = np.zeros((2,4))\n",
    "    ND[0,0] = -0.25*(1-eta) \n",
    "    ND[1,0] = -0.25*(1-xi)\n",
    "    ND[0,1] = 0.25*(1-eta)\n",
    "    ND[1,1] = -0.25*(1+xi)\n",
    "    ND[0,2] = 0.25*(1+eta) \n",
    "    ND[1,2] = 0.25*(1+xi)\n",
    "    ND[0,3] = -0.25*(1+eta) \n",
    "    ND[1,3] = 0.25*(1-xi)\n",
    "    return ND\n",
    "\n",
    "def jacobian(ND,x):\n",
    "    jcb_tr = np.dot(ND,x) # this is transpose(J)\n",
    "    invJ = np.linalg.inv(jcb_tr) # note that invJ = inverse matrix of transpose(J)\n",
    "    detJ = np.linalg.det(jcb_tr) # note that det(J) = det(transpose(J))\n",
    "    return invJ,detJ\n",
    "\n",
    "def local_a(x,ngp,xg,wg):\n",
    "    a = np.zeros((4,4))\n",
    "    for i in range(ngp):\n",
    "        ND = dNdrs(xg[i,0],xg[i,1])\n",
    "        [invJ,detJ] = jacobian(ND,x)\n",
    "        temp1 = np.matmul(invJ,ND)\n",
    "        temp2 = np.transpose(temp1)\n",
    "        a += np.matmul(temp2,temp1)*wg[i]*detJ\n",
    "    return a\n",
    "\n",
    "def assemble(ne,n,ln,ndir,idir,x,ngp,xg,wg):\n",
    "    Keff = np.zeros((n,n))\n",
    "    for i in range(ne):        \n",
    "        # get local matrices of an element\n",
    "        loc_a = local_a(x[ln[i,:],:],ngp,xg,wg)\n",
    "        \n",
    "        # assemble local matrices\n",
    "        for j in range(4):\n",
    "            for k in range(4):\n",
    "                Keff[ln[i][j]][ln[i][k]] += loc_a[j][k]\n",
    "        \n",
    "    # swap elements to apply Dirichlet BC\n",
    "    for i in idir:\n",
    "        Keff[i][:] = 0.\n",
    "        Keff[i][i] = 1.\n",
    "    return Keff\n",
    "\n",
    "def local_cd(x,ngp,xg,wg):\n",
    "    c = np.zeros((4,4))\n",
    "    d1 = np.zeros((4,4))\n",
    "    d2 = np.zeros((4,4))    \n",
    "    for i in range(ngp):\n",
    "        sha = Shapefunctions(xg[i,0],xg[i,1])\n",
    "        ND = dNdrs(xg[i,0],xg[i,1])\n",
    "        [invJ,detJ] = jacobian(ND,x)\n",
    "        temp = np.matmul(invJ,ND)\n",
    "        for j in range(4):\n",
    "            for k in range(4):\n",
    "                c[j][k] += sha[j]*sha[k]*wg[i]*detJ\n",
    "                d1[j][k] += sha[j]*temp[0,k]*wg[i]*detJ\n",
    "                d2[j][k] += sha[j]*temp[1,k]*wg[i]*detJ\n",
    "    return c,d1,d2\n",
    "\n",
    "def v_Galerkin(n,ne,ln,ngp,xg,wg,x,phi):\n",
    "    unod = np.zeros((n,2))\n",
    "    C = np.zeros((n,n))\n",
    "    D1 = np.zeros((n,n))\n",
    "    D2 = np.zeros((n,n))\n",
    "    for i in range(ne):\n",
    "        # get local matrices of an element\n",
    "        [loc_c,loc_d1,loc_d2] = local_cd(x[ln[i,:],:],ngp,xg,wg)\n",
    "        #print(loc_d1)\n",
    "        # assemble local matrices\n",
    "        for j in range(4):\n",
    "            for k in range(4):\n",
    "                C[ln[i][j]][ln[i][k]] += loc_c[j][k]\n",
    "                D1[ln[i][j]][ln[i][k]] += loc_d1[j][k]\n",
    "                D2[ln[i][j]][ln[i][k]] += loc_d2[j][k]\n",
    "                \n",
    "    # solve u\n",
    "    b = np.matmul(D1,phi)\n",
    "    unod[:,0] = np.linalg.solve(C,b)\n",
    "    # solve v\n",
    "    b = np.matmul(D2,phi)\n",
    "    unod[:,1] = np.linalg.solve(C,b)\n",
    "    return unod\n",
    "\n",
    "def v_average(n,ne,ln,xc,x,sol):   \n",
    "    unod = np.zeros((n,2))\n",
    "    c = np.zeros(n)    \n",
    "    phi = np.zeros(4)\n",
    "    for i in range(ne):\n",
    "        c[ln[i][0]] += 1.\n",
    "        c[ln[i][1]] += 1.\n",
    "        c[ln[i][2]] += 1.\n",
    "        c[ln[i][3]] += 1.\n",
    "        for j in range(4):\n",
    "            # get inverse of JT and determinant of J\n",
    "            ND = dNdrs(xc[j,0],xc[j,1])\n",
    "            [invJ,detJ] = jacobian(ND,x[ln[i][:]][:])\n",
    "            phi = [sol[ln[i][0]], sol[ln[i][1]], sol[ln[i][2]], sol[ln[i][3]]]\n",
    "            unod[ln[i,j],:] += np.dot(invJ,np.dot(ND,phi))\n",
    "    for i in range(n):\n",
    "        unod[i][0] = unod[i][0]/c[i]\n",
    "        unod[i][1] = unod[i][1]/c[i]\n",
    "    return unod\n",
    "\n",
    "def bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi):\n",
    "    # coordinate transformation and assemble local matrices\n",
    "    Keff = assemble(ne,n,ln,ndir,idir,x,ngp,xg,wg)\n",
    "        \n",
    "    # build load vector that satisfies free-surface BC (Dirichlet phi is given)\n",
    "    F = load_phi(n,ndir,idir,phi)\n",
    "        \n",
    "    # solve potentialin the non-inertia frame\n",
    "    phi = np.linalg.solve(Keff,F)\n",
    "        \n",
    "    # calculate nodal velocity by averaging cell velocity\n",
    "    if mtd_unod == 'avg':\n",
    "        unod = v_average(n,ne,ln,xc,x,phi)\n",
    "    elif mtd_unod == 'glk':\n",
    "        unod = v_Galerkin(n,ne,ln,ngp,xg,wg,x,phi)\n",
    "    \n",
    "    # calculate dpdt and phit on free surface\n",
    "    [phit,dpdt] = FSBC_phit(grav,h,ifs,n,x,unod,acc)\n",
    "        \n",
    "    # return current velocity and dpdt (if ale is used, return mesh velocity and dpdt)\n",
    "    u_cur = unod\n",
    "    dpdt_cur = dpdt\n",
    "        \n",
    "    # calculate mesh velocity, convective velocity, and time derivative of phi in ALE frame\n",
    "    if mtd_remesh == 'ale':\n",
    "        normy = normal_vector(nfs,ifs,n,x)\n",
    "        [vm,cv] = v_mesh(nfs,ifs,normy,unod,0.)\n",
    "        dpdtm = dpdt_mesh(nfs,ifs,dpdt,cv,unod)\n",
    "        for i in range(nfs):\n",
    "            u_cur[ifs[i],:] = vm[i,:]\n",
    "            dpdt_cur[ifs[i]] = dpdtm[i]\n",
    "            \n",
    "    return u_cur,dpdt_cur\n",
    "\n",
    "def initial_FSBC_slo(n,ndir,idir,L,h,x,u0,a0):\n",
    "    phi = np.zeros(n)\n",
    "    phit = np.zeros(n)\n",
    "    for i in idir:\n",
    "        phi[i] = -x[i][0]*u0\n",
    "    return phi,phit\n",
    "\n",
    "def initial_FSBC_std(a,n,ndir,idir,L,h,x):\n",
    "    phi = np.zeros(n)\n",
    "    phit = np.zeros(n)\n",
    "    for i in idir:\n",
    "        phi[i] = a*np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*h/L)\n",
    "    return phi,phit\n",
    "\n",
    "def normal_vector(nfs,ifs,n,x):\n",
    "    le = np.zeros(nfs-1)\n",
    "    ne = np.zeros((nfs-1,2))\n",
    "    normy = np.zeros((nfs,2))\n",
    "    \n",
    "    # normal vector of FS side of element\n",
    "    for i in range(nfs-1): \n",
    "        dx = x[ifs[i+1]][0]-x[ifs[i]][0]\n",
    "        dy = x[ifs[i+1]][1]-x[ifs[i]][1]\n",
    "        le[i] = np.sqrt(dx**2+dy**2)\n",
    "        ne[i][0] = -dy/le[i]\n",
    "        ne[i][1] = dx/le[i]\n",
    "    \n",
    "    # nodal normal vector\n",
    "    normy[0][0] = ne[0][0]\n",
    "    normy[0][1] = ne[0][1]\n",
    "    normy[-1][0] = ne[-1][0]\n",
    "    normy[-1][1] = ne[-1][1]\n",
    "    for i in range(1,nfs-1): \n",
    "        normy[i][0] = (le[i-1]*ne[i-1][0]+le[i]*ne[i][0])/(le[i-1]+le[i])\n",
    "        normy[i][1] = (le[i-1]*ne[i-1][1]+le[i]*ne[i][1])/(le[i-1]+le[i])         \n",
    "    return normy\n",
    "\n",
    "def v_mesh(nfs,ifs,normy,unod,vt):\n",
    "    vm = np.zeros((nfs,2))\n",
    "    cv = np.zeros((nfs,2))\n",
    "    for i in range(nfs):\n",
    "        vm[i][0] = 0. # the mesh velocity has no horizontal velocity\n",
    "        vm[i][1] = unod[ifs[i]][1]+(unod[ifs[i]][0]-vm[i][0])*normy[i][0]/normy[i][1]\n",
    "        cv[i][0] = unod[ifs[i]][0]-vm[i][0]\n",
    "        cv[i][1] = unod[ifs[i]][1]-vm[i][1]\n",
    "    return vm,cv\n",
    "\n",
    "def dpdt_mesh(nfs,ifs,dpdt,cv,unod):\n",
    "    dpdtm = np.zeros(nfs)\n",
    "    for i in range(nfs):\n",
    "        dpdtm[i] = dpdt[ifs[i]]-cv[i][0]*unod[ifs[i]][0]-cv[i][1]*unod[ifs[i]][1]\n",
    "    return dpdtm\n",
    "\n",
    "def FSBC_phit(grav,h,ifs,n,x,unod,acc):\n",
    "    dpdt = np.zeros(n)\n",
    "    phit = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        dpdt[i] = 0.5*(unod[i][0]**2+unod[i][1]**2)-grav*(x[i][1]-h)-x[i][0]*acc\n",
    "        phit[i] = dpdt[i]-(unod[i][0]**2+unod[i][1]**2)\n",
    "    return phit,dpdt\n",
    "\n",
    "def total_phi(n,x,phi,u,v):\n",
    "    PHI = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        PHI[i] = phi[i] + x[i][0]*u + x[i][1]*v\n",
    "    return PHI\n",
    "\n",
    "def moving_phi(n,x,PHI,u,v):\n",
    "    phi = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        phi[i] = PHI[i] - x[i][0]*u - x[i][1]*v\n",
    "    return phi\n",
    "\n",
    "def remesh_EL(mtd_mesh,dt,nx,ny,n,ifs,x,unod,phi,dpdt):\n",
    "    # update free-surface node and potential\n",
    "    new_x = np.zeros((n,2))\n",
    "    new_phi = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        new_x[i][0] = x[i][0]+dt*unod[i][0]\n",
    "        new_x[i][1] = x[i][1]+dt*unod[i][1]\n",
    "        new_phi[i] = phi[i]+dt*dpdt[i]\n",
    "        \n",
    "    # rezoneing all nodes\n",
    "    for i in range(nx):\n",
    "        if mtd_mesh == 'even':\n",
    "            yi = np.linspace(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='cos':\n",
    "            yi = x_cos(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='pwr':\n",
    "            yi = x_power(0.,new_x[ifs[i]][1],ny)\n",
    "        for j in range(ny):\n",
    "            c = j*nx+i\n",
    "            new_x[c][0] = new_x[ifs[i]][0]\n",
    "            new_x[c][1] = yi[j]\n",
    "    return new_x,new_phi\n",
    "\n",
    "def remesh_EL_AB2(mtd_mesh,dt,nx,ny,n,ifs,x,unod,unod2,phi,dpdt,dpdt2):\n",
    "    # unod is the velocity at kth time step\n",
    "    # unod2 is the velocity at (k-1)th time step\n",
    "    # update the position and potential by AB2 method\n",
    "    new_x = np.zeros((n,2))\n",
    "    new_phi = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        new_x[i][0] = x[i][0]+0.5*dt*(3.*unod[i][0]-unod2[i][0])\n",
    "        new_x[i][1] = x[i][1]+0.5*dt*(3.*unod[i][1]-unod2[i][1])\n",
    "        new_phi[i] = phi[i]+0.5*dt*(3.*dpdt[i]-dpdt2[i])\n",
    "        \n",
    "    # rezoneing all nodes\n",
    "    for i in range(nx):\n",
    "        if mtd_mesh == 'even':\n",
    "            yi = np.linspace(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='cos':\n",
    "            yi = x_cos(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='pwr':\n",
    "            yi = x_power(0.,new_x[ifs[i]][1],ny)\n",
    "        for j in range(ny):\n",
    "            c = j*nx+i\n",
    "            new_x[c][0] = new_x[ifs[i]][0]\n",
    "            new_x[c][1] = yi[j]\n",
    "    return new_x,new_phi\n",
    "\n",
    "def remesh_ALE(mtd_mesh,dt,nx,ny,n,nfs,ifs,x,vm,phi,dpdtm):\n",
    "    # update free-surface node and potential\n",
    "    new_x = np.zeros((n,2))\n",
    "    new_phi = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        new_x[i,0] = x[i,0]+dt*vm[i][0]\n",
    "        new_x[i,1] = x[i,1]+dt*vm[i][1]\n",
    "        new_phi[i] = phi[i]+dt*dpdtm[i]\n",
    "        \n",
    "    # rezoneing all nodes\n",
    "    for i in range(nx):\n",
    "        if mtd_mesh == 'even':\n",
    "            yi = np.linspace(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='cos':\n",
    "            yi = x_cos(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='pwr':\n",
    "            yi = x_power(0.,new_x[ifs[i]][1],ny)\n",
    "        for j in range(ny):\n",
    "            c = j*nx+i\n",
    "            new_x[c][0] = new_x[ifs[i]][0]\n",
    "            new_x[c][1] = yi[j]\n",
    "\n",
    "    return new_x,new_phi\n",
    "\n",
    "def remesh_ALE_AB2(mtd_mesh,dt,nx,ny,n,nfs,ifs,x,vm,vm2,phi,dpdtm,dpdtm2):\n",
    "    # update free-surface node and potential\n",
    "    new_x = np.zeros((n,2))\n",
    "    new_phi = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        new_x[i,0] = x[i,0]+0.5*dt*(3.*vm[i][0]-vm2[i][0])\n",
    "        new_x[i,1] = x[i,1]+0.5*dt*(3.*vm[i][1]-vm2[i][1])\n",
    "        new_phi[i] = phi[i]+0.5*dt*(3.*dpdtm[i]-dpdtm2[i])\n",
    "\n",
    "    # rezoneing all nodes\n",
    "    for i in range(nx):\n",
    "        if mtd_mesh == 'even':\n",
    "            yi = np.linspace(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='cos':\n",
    "            yi = x_cos(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='pwr':\n",
    "            yi = x_power(0.,new_x[ifs[i]][1],ny)\n",
    "        for j in range(ny):\n",
    "            c = j*nx+i\n",
    "            new_x[c][0] = new_x[ifs[i]][0]\n",
    "            new_x[c][1] = yi[j]\n",
    "            \n",
    "    return new_x,new_phi\n",
    "\n",
    "def rk4(smth_fs,mtd_mesh,mtd_unod,mtd_remesh,dt,acc0,accn,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi):\n",
    "    # define acc for intermediate step, acc0, acci, accn\n",
    "    acci = 0.5*(acc0+accn)\n",
    "    \n",
    "    if mtd_remesh == 'el': # Lagrangian remesh\n",
    "        # ki = nodal velocity or dpdt, get k1\n",
    "        [u1,dpdt1] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc0,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)\n",
    "        [x1,phi1] = remesh_EL(mtd_mesh,0.5*dt,nx,ny,n,ifs,x,u1,phi,dpdt1)\n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k2\n",
    "        [u2,dpdt2] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acci,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x1,phi1)    \n",
    "        [x2,phi2] = remesh_EL(mtd_mesh,0.5*dt,nx,ny,n,ifs,x,u2,phi,dpdt2)\n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k3\n",
    "        [u3,dpdt3] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acci,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x2,phi2)\n",
    "        [x3,phi3] = remesh_EL(mtd_mesh,dt,nx,ny,n,ifs,x,u3,phi,dpdt3)    \n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k4\n",
    "        [u4,dpdt4] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,accn,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x3,phi3)    \n",
    "        \n",
    "    if mtd_remesh == 'ale': # ALE remesh\n",
    "        # ki = nodal velocity or dpdt, get k1\n",
    "        [u1,dpdt1] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc0,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)         \n",
    "        [x1,phi1] = remesh_ALE(mtd_mesh,0.5*dt,nx,ny,n,nfs,ifs,x,u1,phi,dpdt1)\n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k2\n",
    "        [u2,dpdt2] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acci,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x1,phi1)    \n",
    "        [x2,phi2] = remesh_ALE(mtd_mesh,0.5*dt,nx,ny,n,nfs,ifs,x,u2,phi,dpdt2)\n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k3\n",
    "        [u3,dpdt3] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acci,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x2,phi2)\n",
    "        [x3,phi3] = remesh_ALE(mtd_mesh,dt,nx,ny,n,nfs,ifs,x,u3,phi,dpdt3)    \n",
    "    \n",
    "        # ki = nodal velocity or dpdt, get k4\n",
    "        [u4,dpdt4] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,accn,grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x3,phi3)    \n",
    "    \n",
    "    # define ki for x and phi on the free surface\n",
    "    k1x = []; k1y = []; k1p = []; k2x = []; k2y = []; k2p = []\n",
    "    k3x = []; k3y = []; k3p = []; k4x = []; k4y = []; k4p = []\n",
    "    for i in range(nfs):\n",
    "        k1x.append(u1[ifs[i],0])\n",
    "        k1y.append(u1[ifs[i],1])\n",
    "        k1p.append(dpdt1[ifs[i]])\n",
    "        k2x.append(u2[ifs[i],0])\n",
    "        k2y.append(u2[ifs[i],1])\n",
    "        k2p.append(dpdt2[ifs[i]])        \n",
    "        k3x.append(u3[ifs[i],0])\n",
    "        k3y.append(u3[ifs[i],1])\n",
    "        k3p.append(dpdt3[ifs[i]])\n",
    "        k4x.append(u4[ifs[i],0])\n",
    "        k4y.append(u4[ifs[i],1])\n",
    "        k4p.append(dpdt4[ifs[i]])  \n",
    "\n",
    "    [new_x,new_phi] = remesh_rk4(smth_fs,mtd_mesh,dt,nx,ny,n,nfs,ifs,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y,k1p,k2p,k3p,k4p,x,phi)\n",
    "    return new_x,new_phi\n",
    "\n",
    "def remesh_rk4(smth_fs,mtd_mesh,dt,nx,ny,n,nfs,ifs,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y,k1p,k2p,k3p,k4p,x,phi):  \n",
    "    # update free-surface node and potential\n",
    "    new_x = np.zeros((n,2))\n",
    "    new_phi = np.zeros(n)\n",
    "    for i in range(nfs):\n",
    "        new_x[ifs[i]][0] = x[ifs[i]][0]+dt*(k1x[i]+2.*k2x[i]+2.*k3x[i]+k4x[i])/6.\n",
    "        new_x[ifs[i]][1] = x[ifs[i]][1]+dt*(k1y[i]+2.*k2y[i]+2.*k3y[i]+k4y[i])/6.\n",
    "        new_phi[ifs[i]] = phi[ifs[i]]+dt*(k1p[i]+2.*k2p[i]+2.*k3p[i]+k4p[i])/6.\n",
    "\n",
    "    if smth_fs:\n",
    "        [new_x,new_phi] = smooth(nfs,ifs,new_x,new_phi)\n",
    "        \n",
    "    # rezoneing all nodes\n",
    "    for i in range(nx):\n",
    "        if mtd_mesh == 'even':\n",
    "            yi = np.linspace(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='cos':\n",
    "            yi = x_cos(0.,new_x[ifs[i]][1],ny)\n",
    "        elif mtd_mesh =='pwr':\n",
    "            yi = x_power(0.,new_x[ifs[i]][1],ny)\n",
    "        for j in range(ny):\n",
    "            c = j*nx+i\n",
    "            new_x[c][0] = new_x[ifs[i]][0]\n",
    "            new_x[c][1] = yi[j]\n",
    "    return new_x,new_phi\n",
    "\n",
    "def pressure(tho,grav,h,ic,ifs,n,x,dpdt_cur,u_cur):\n",
    "    phit = np.zeros(n)\n",
    "    for i in ifs:\n",
    "        phit[i] = dpdt_cur[i]-(u_cur[i,0]**2+u_cur[i,1]**2)\n",
    "    \n",
    "    # find the atmospheric pressure (pressure on the free surface)\n",
    "    patm = tho*(phit[ifs[ic]]+0.5*(u_cur[ifs[ic]][0]**2+u_cur[ifs[ic]][1]**2)+grav*(x[ifs[ic],1]-h))\n",
    "    p1 = np.zeros(n)\n",
    "    p2 = np.zeros(n)\n",
    "    p3 = np.zeros(n)\n",
    "    p = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        p1[i] = tho*phit[i]\n",
    "        p2[i] = tho*0.5*(u_cur[i][0]**2+u_cur[i][1]**2)\n",
    "        p3[i] = tho*grav*(x[i][1]-h)\n",
    "        p[i] = p1[i]+p2[i]+p3[i]-patm\n",
    "    return patm,p1,p2,p3,p\n",
    "\n",
    "def fs_animation(t,x,y,L,h):\n",
    "    display.clear_output(wait=True)\n",
    "    plt.plot(x,y,'-ob')\n",
    "    plt.xlim([-0.05, L+0.05])\n",
    "    plt.ylim([h-0.05, h+0.05])\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('y')\n",
    "    plt.title(str(t))\n",
    "    plt.show()\n",
    "\n",
    "def arc_length(nfs,ifs,x,phi):\n",
    "    # raw free-surface data\n",
    "    fx = []\n",
    "    fy = []\n",
    "    for i in ifs:\n",
    "        fx.append(x[i][0])\n",
    "        fy.append(x[i][1])\n",
    "    \n",
    "    # calculate arc length\n",
    "    cs = interp.CubicSpline(fx,fy)\n",
    "    csder = cs.derivative(1)\n",
    "    ns = 5*nfs\n",
    "    xs = np.linspace(fx[0],fx[-1],ns)\n",
    "    dx = xs[1]-xs[0]\n",
    "    s = 0.5*dx*(np.sqrt(1.+csder(xs[0])**2)+np.sqrt(1.+csder(xs[-1])**2))\n",
    "    for i in range(ns-2):\n",
    "        s += dx*np.sqrt(1.+csder(xs[i+1])**2)\n",
    "    return s\n",
    "\n",
    "def cal_mass(ne,ln,x,area_exa):\n",
    "    a = np.zeros(2)\n",
    "    b = np.zeros(2)\n",
    "    c = np.zeros(2)\n",
    "    d = np.zeros(2)\n",
    "    area = 0.\n",
    "    for i in range(ne):\n",
    "        for j in range(2):\n",
    "            a[j] = x[ln[i,1],j]-x[ln[i,0],j]\n",
    "            b[j] = x[ln[i,3],j]-x[ln[i,0],j]\n",
    "            c[j] = x[ln[i,1],j]-x[ln[i,2],j]\n",
    "            d[j] = x[ln[i,3],j]-x[ln[i,2],j]\n",
    "        da1 = 0.5*np.abs(np.cross(a,b))\n",
    "        da2 = 0.5*np.abs(np.cross(c,d))\n",
    "        area += da1+da2\n",
    "    err = (area-area_exa)/area_exa*100.\n",
    "    return err\n",
    "\n",
    "def guass_qua(gq_order):\n",
    "    if gq_order == 1:\n",
    "        ng = 1\n",
    "    elif gq_order == 2:\n",
    "        ng = 2\n",
    "    elif gq_order == 3:\n",
    "        ng = 3\n",
    "    elif gq_order == 4:\n",
    "        ng = 4\n",
    "    \n",
    "    # distribute ng-by-ng node on the local element\n",
    "    [x,w] = special.p_roots(ng)\n",
    "    ngp = ng*ng\n",
    "    xg = np.zeros((ngp,2))\n",
    "    wg = []\n",
    "    c = 0\n",
    "    for i in range(ng):\n",
    "        for j in range(ng):\n",
    "            xg[c][0] = x[j]\n",
    "            xg[c][1] = x[i]\n",
    "            wg.append(w[i]*w[j])\n",
    "            c += 1\n",
    "    return ngp,xg,wg\n",
    "\n",
    "def main(x,phi,phit):\n",
    "    # TIME MARCHING START\n",
    "    c_frame = 0\n",
    "    for i in range(nt-1):\n",
    "#        fs_animation(time[i],x[ifs[:],0],x[ifs[:],1],L,h)\n",
    "\n",
    "        # output wave elevation on both side walls\n",
    "        ita1 = x[ifs[0],1]-h\n",
    "        ita2 = x[ifs[-1],1]-h\n",
    "        print(time[i],ita1,ita2,file=wav)\n",
    "        \n",
    "        # check to avoid high wave elevation\n",
    "        check_ita(ita1,h)\n",
    "        \n",
    "        # calculate the numerical mass err in %\n",
    "        m_err = cal_mass(ne,ln,x,Vw)\n",
    "        print(time[i],m_err,file=mas)\n",
    "        \n",
    "        # current velocity and pressure distribution\n",
    "        [u_cur,dpdt_cur] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc[i],grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)\n",
    "        [patm,p1,p2,p3,p] = pressure(tho,grav,h,ic,ifs,n,x,dpdt_cur,u_cur)\n",
    "        if i % n_frame == 0:\n",
    "            save_pic(c_frame,ln,x,u_cur,p)\n",
    "            c_frame += 1\n",
    "        \n",
    "        # time integration\n",
    "        if mtd_intgr == 'fir': # 1st-order explicit\n",
    "            [u_cur,dpdt_cur] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc[i],grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)\n",
    "            if mtd_remesh == 'el': # Lagrangian remesh\n",
    "                [new_x,new_phi] = remesh_EL(mtd_mesh,dt,nx,ny,n,ifs,x,u_cur,phi,dpdt_cur)\n",
    "            elif mtd_remesh == 'ale': # ALE remesh\n",
    "                [new_x,new_phi] = remesh_ALE(mtd_mesh,dt,nx,ny,n,nfs,ifs,x,u_cur,phi,dpdt_cur)            \n",
    "            \n",
    "        elif mtd_intgr == 'ab2': # AB2\n",
    "            [u_cur,dpdt_cur] = bvp_solver(mtd_mesh,mtd_unod,mtd_remesh,acc[i],grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)     \n",
    "            if mtd_remesh == 'el': # Lagrangian remesh with AB2\n",
    "                [new_x,new_phi] = remesh_EL_AB2(mtd_mesh,dt,nx,ny,n,ifs,x,u_cur,u_pre,phi,dpdt_cur,dpdt_pre)\n",
    "            elif mtd_remesh == 'ale': # ALE remesh with AB2\n",
    "                [new_x,new_phi] = remesh_ALE_AB2(mtd_mesh,dt,nx,ny,n,nfs,ifs,x,u_cur,u_pre,phi,dpdt_cur,dpdt_pre)\n",
    "            u_pre = u_cur\n",
    "            dpdt_pre = dpdt_cur\n",
    "            \n",
    "        elif mtd_intgr == 'rk4': # rk4\n",
    "            [new_x,new_phi] = rk4(smth_fs,mtd_mesh,mtd_unod,mtd_remesh,dt,acc[i],acc[i+1],grav,h,ne,nx,ny,n,ln,ndir,idir,nfs,ifs,ngp,xg,wg,xc,x,phi)\n",
    "\n",
    "        # update the free surface\n",
    "        x = new_x\n",
    "        phi = new_phi\n",
    "        \n",
    "        print(i)\n",
    "\n",
    "def save_pic(step,ln,x,u_cur,p):\n",
    "    plt.figure(figsize=(16,4)) #; plt.rcParams.update({'font.size': 22})\n",
    "    for eN, e in enumerate(ln):\n",
    "        element_nodes = x[e[0]], x[e[1]], x[e[2]], x[e[3]]\n",
    "        X = [n[0] for n in element_nodes]; Y = [n[1] for n in element_nodes]\n",
    "        X.append(element_nodes[0][0]); Y.append(element_nodes[0][1])\n",
    "        plt.plot(X,Y,'-k',linewidth=0.5) #;plt.savefig('C{0:04d}.png'.format(step))\n",
    "    \n",
    "    plt.tricontourf(x[:,0], x[:,1], p, 20)\n",
    "    plt.colorbar()\n",
    "    plt.quiver(x[:,0],x[:,1],u_cur[:,0],u_cur[:,1],scale=1,width=0.002)\n",
    "    #plt.show()\n",
    "    plt.savefig('C{0:04d}.png'.format(step))\n",
    "    \n",
    "#************************************************************************\n",
    "#   The ALE-FEM for 2D sloshing problem in potential flow\n",
    "#   Wen-Huai Tsao 2022 Summer at LSU\n",
    "#************************************************************************\n",
    "import numpy as np\n",
    "import math\n",
    "import random as rd\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "from scipy.sparse import csr_matrix,linalg\n",
    "import scipy.interpolate as interp\n",
    "from scipy import special\n",
    "from scipy import integrate\n",
    "#from scipy.sparse.linalg import spsolve\n",
    "#from cython.parallel import prange\n",
    "from IPython import display\n",
    "\n",
    "# OPEN FILES\n",
    "#************************************************************************\n",
    "wav = open(\"wav_quad.dat\", \"w\")\n",
    "mas = open(\"mass_quad.dat\", \"w\")\n",
    "test = open(\"test.txt\", \"w\")\n",
    "\n",
    "# PHYSICS CONFIGURATIONS\n",
    "#************************************************************************\n",
    "# physical constant\n",
    "tho = 998.2\n",
    "grav = 9.81\n",
    "\n",
    "# tank dimensions and 1st-order natural frequency (rad/s)\n",
    "L = 0.8\n",
    "h = 0.1\n",
    "Vw = L*h\n",
    "wn = math.sqrt(np.pi*grav*np.tanh(np.pi*h/L)/L)\n",
    "\n",
    "# claim for some variables to be saved\n",
    "#u_cur = np.zeros((n,2))\n",
    "#dpdt_cur = np.zeros(n)\n",
    "#u_pre = np.zeros((n,2))\n",
    "#dpdt_pre = np.zeros(n)\n",
    "\n",
    "# NUMERICAL CONFIGURATIONS\n",
    "#************************************************************************\n",
    "# number of threads\n",
    "nthrd = 8\n",
    "    \n",
    "# n-by-n gaussian quadrature nodes on a element\n",
    "gq_order = 2 #  n = 1, 2, 3, 4\n",
    "\n",
    "# Guassian quadrature\n",
    "[ngp,xg,wg] = guass_qua(gq_order)    \n",
    "\n",
    "# for nodal velocity by average method\n",
    "xc = loc_nod()\n",
    "\n",
    "# mesh size and time interval\n",
    "he = 0.01\n",
    "dt = 0.005\n",
    "T = 20.\n",
    "nt = int(T/dt)+1\n",
    "time = np.linspace(0.,T,nt)\n",
    "\n",
    "# ifs[ic] is the index of free-surce node for Bernoulli constant\n",
    "ic = int(L/he/2)+1\n",
    "\n",
    "# save contour and quiver every n frames\n",
    "n_frame = 20\n",
    "\n",
    "# INITIAL MESH\n",
    "#************************************************************************\n",
    "# node distribution along vertical direction\n",
    "mtd_mesh = 'pwr' # 'even' or 'cos' or 'pwr'\n",
    "\n",
    "# mesh and l2g index\n",
    "[nex,ney,ne,nx,ny,n,ln,ndir,idir,nneu1,ineu1,nneu2,ineu2,nneu3,ineu3,ifs,nfs,x] = mesh(L,h,he,mtd_mesh)\n",
    "#plotMesh(ln,x)\n",
    "\n",
    "# SELECT METHODS\n",
    "#************************************************************************\n",
    "mtd_unod = 'glk' # 'avg' or 'glk'\n",
    "mtd_remesh = 'ale' # 'el' or 'ale'\n",
    "mtd_intgr = 'rk4' # 'fir' or 'ab2' or 'rk4'\n",
    "smth_fs = False # True or False\n",
    "\n",
    "# SELECT CASES\n",
    "#************************************************************************\n",
    "# std = standing wave in stationary tank; slo = sloshing in moving tank\n",
    "exp_name = 'slo'\n",
    "\n",
    "# amplitude of potential on free surface\n",
    "ap = 0.01\n",
    "\n",
    "# ground excitation = af*sin(wf*t)\n",
    "beta = 1.0\n",
    "wf = beta*wn\n",
    "af = 0.0005\n",
    "\n",
    "# define tank motion\n",
    "dis = np.zeros(nt)\n",
    "vel = np.zeros(nt)\n",
    "acc = np.zeros(nt)\n",
    "if exp_name == 'std': # standing wave by giving initial velocity on free surface\n",
    "    [phi,phit] = initial_FSBC_std(ap,n,ndir,idir,L,h,x)\n",
    "elif exp_name == 'slo': # harmonic excitation: x = af*sin(wf*t)\n",
    "    for i in range(nt):\n",
    "        dis[i] = af*np.sin(wf*dt*i)\n",
    "        vel[i] = wf*af*np.cos(wf*dt*i)\n",
    "        acc[i] = -wf*wf*af*np.sin(wf*dt*i)\n",
    "    [phi,phit] = initial_FSBC_slo(n,ndir,idir,L,h,x,vel[0],vel[0])\n",
    "\n",
    "# RUN SIMULATION\n",
    "#************************************************************************\n",
    "%time main(x,phi,phit) # Q: why x is not a global variable\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8484d1fe-649d-48d7-a15b-a92f5a03affe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'mass error (%)')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEGCAYAAADIRPqpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABSEElEQVR4nO29d3gc53Xo/TuLjkXvhSAJkiApNlESRfXe3aQkdmwltmTHJe52cpMb5UuuldhfEvtznOJcx4msyJaur+0osSzJtmSZ6pZESqREik3sBSB6X9RF2ff7Y2dAEFxgC2Z2luD5Pc8+mPLOO2dIYM6e854ixhgURVEUxUt8XgugKIqiKKqMFEVRFM9RZaQoiqJ4jiojRVEUxXNUGSmKoiiek+61AOcqZWVlZunSpV6LoSiKcs7w5ptvdhljyiOdU2WUIEuXLmXHjh1ei6EoinLOICInZzunbjpFURTFc1QZKYqiKJ6jykhRFEXxHE+VkYjcLiIHReSIiNwX4XyWiPyndf51EVk67dyfW8cPisht0eYUkXprjiPWnJnR7qEoiqIkB8+UkYikAd8B7gDWAHeLyJoZwz4O9BpjVgD/CHzDunYN8CFgLXA78K8ikhZlzm8A/2jN1WvNPes9FEVRlOThpWW0GThijDlmjBkDfgLcOWPMncDD1vZ/AzeJiFjHf2KMCRpjjgNHrPkizmldc6M1B9acd0W5h6IoipIkvFRGtUDTtP1T1rGIY4wxE0A/UDrHtbMdLwX6rDlm3mu2e5yFiHxKRHaIyI7Ozs6YH1RRFEWZG80zigNjzAPAAwCbNm0653pvNPUM8+v97YxPhrh0aQmXLCn2WqSYMcbwVmMfI2OTbK4vITNdY28UZSHhpTJqBuqm7S+yjkUac0pE0oFCoDvKtZGOdwNFIpJuWT/Tx892jwXFf+1o4v/52R7GJ0/r0Ls21vDND1xIRlpqv9hHxyf5wo93smV/OwArK/P4tw9fwrLyPI8li403jvfwN0+9w9hEiC/f3MBta6u8Fiku2gOjDAYnWFbmRz3Yilt4+RbaDjRYUW6ZhAMSnpwx5kngXmv7/cDzJtwN8EngQ1YkXD3QALwx25zWNS9Yc2DN+USUeywYnj/Qzv/86W4uqy/llT+7gV1fuYUv3dTA47ta+LunDngtXlT++uf72LK/nfvuWM2/3H0RXYNj3PPQG/QOjXktWlQOtAW456HX6RkKMhkK8ZkfvsmLBzu8FitmHtl6giv+7jlu+tZL/I9H3yYUOrf+NDoGRjnSMeC1GEoMeKaMLAvl88AzwDvAo8aYfSLyVRF5nzXsP4BSETkC/DFwn3XtPuBRYD/wK+BzxpjJ2ea05voz4I+tuUqtuWe9x0KhpW+EL/14F2uqC/jePZtYVJxLUW4mf3TLSu69YgkPvXqcN473eC3mrLx0qJMfv9HEZ65fzqevW857L6zhoY9eSmv/KN/actBr8ebEGMNXnthHbmY6P/3MlTz+uatoqMjnzx/bw+j4pNfiRWVvcz9/9eQ+rl9VwSeuruexnc388PVZq7mkHM/ub+eab7zAzf/wMn/31DteixMXxhh+va+NR7c3nXNfABJFFpgRkDQ2bdpkUr02nTGGj/1gO68f6+GZL1/L4tLcM84Pj01w87deoiAng1984WrSU8xdNz4Z4vZ/epnJkOGZP7qWrPS0qXP3P7GX/7PtJC/+yQ1nPVeq8OLBDj76/e187a51fOTyJQC8dqSL33vwdf72t9bze5ct9ljCufnMD9/k1SNdvHLfjeRnpXP397ZxrHOI3/zZDWf8X6Qi/cPjXP3/Pc+S0lzqy/L4+dstPPqHV7C5vsRr0WLika0n+MoT4e/RH71yKX/1vrUeS+QMIvKmMWZTpHOp9fZRHOW5dzp48WAnf3rbqogv7NzMdP7i3Ws40DbAL3a3eiDh3Dy5q4WjnUP8xbvXnPXy+8z1K/CJ8H+2nfBGuBj4wWsnqMjP4oObTi9jXrG8lA2LCnnwlWOk8hfBU73DPLOvjQ9fvoSC7AxEhE9ft5yOgSAvHEj9SNIfvHaCgdEJvvn+C/nm+zdQ6s/kgZePei1WTIyMTfLNZw5yTUMZd2+u45GtJzjaOei1WK6jymiBYozhH7YcYnFJLh+5Ysms4+5YV8Wqynz+9wtHUsodYIzh+68dZ0VFHjdfUHHW+arCbG5bV8WjO04RnEg9l1dTzzAvHuzk7s2Lz4j8ExE+fNkSjnUO8fapfg8lnJun97QRMnD35tPW29Uryij1Z/Lz3S0eShadUMjww9dPcuPqCi6oLiA7I40PXlrH8wc66DkH1hmf2tPKwOgEn79hBX98yyp8IvzkjUavxXIdVUYLlOcPdLC/NcAXb2qYM1rO5xM+d+MKjnQM8vyB1FlYf6uxl73NAT565dJZI7g+cMki+kfGeelg6n1Tf2pP2NJ8/yWLzjp3+/oqMtN9PL5zZvBo6vCrfW2srSmgruS0RZ2e5uNd66t57p32lF7zequxl86BIHdurJk6dse6akIGnnun3UPJYuOxnadYWprL5voSyvOzuG5lOU/vbUtpS9oJVBktUB7ZepLKgqwz/iBn4451VZTlZfGfO5qijk0Wj73VTE5GGr910cw86NNcZX1Tf+Lt1Pum/ss9rWxYVHjGy9ymIDuD61aWs2V/e0q+YDoCo7x5spfbI4Sg37i6gtHxEDtO9HogWWz8am8bmWk+blx92qJeV1tAbVEOv96f2spoZGyS7cd7uWVN5dSXsOtXlXOqd4TjXUMeS+cuqowWICe7h3jpUNhFFEsOUUaaj9+5pJbnD3TQMTCaBAnnZmIyxNN727jpggr8WbOnwmWk+Xj3hmqe3d/OUHBi1nHJpqlnmN2n+nnX+upZx9ywqoLmvhGOdKTeWsALVuj5rRGU0WXLSshIE35zJPWsUZtXjnRxaX0x+dkZU8dEhGtXlrHtWHdKuaNn8vrxbsYmQ1zdcLoZ6nUrw0r15UOp+2/uBKqMFiA/er2RNJ+c4e+Pxu9uqmMyZHhyl/dWxmtHu+kZGuM9G2Kx6qoJToT4zeGuJEgWG3Ye0VzJrdetCr9sXkrBF8y2Yz2U5WWysvLspOLczHQuXlzMKyn07z2d/pFxDrYPsHnp2RW9Ll1awsDoBAfbUzfv6LWj3WSm+di89HTU3+LSXBYV57A9ha1RJ1BltMAYHZ/kP3c0cdvaSioLsmO+bnl5HmtrCqbWOrzkF7tbyMtK5/pV5VHHblpaTEF2ekqtBbxypItFxTksnSPkvLYoh4aKPF5MsfUuYwzbjnVz2bLSWdfqLqsv4Z3WAMNjqWON2rx5sgdjiBjCfan1gt9xInXz6nY19bG2toCczDOjRy+sK+LtU33eCJUkVBktMH65u5W+4XE+fPnsEXSz8a711bzV2Edr/4gLksXG2ESIX+1t49Y1lWRnRM9lyUjzcd2qCl442JES7peJyRCvHe3m6hVlUUvnXNNQzvYTPYxNhJIkXXQae4Zp7R/l8mURawUDcNHiYkIG3m5KvWjA7Sd6SfcJG+uKzjq3qDiHivws3jyZmhZGKGTY3xJgfW3hWec21BZyqnfknIgGTBRVRguM/3qzifoyP1fM8TKZjTvWhd1Kv9rb5rRYMbPjRA+B0QluXxd7/babL6iga3CM3c3evxz3NPczMDrB1Q1lUcdeurSY4ESIvS3ey22z7Vi4LOMVy2ZPDrVf9DubUu+lvre5n5WV+WdZFhBeN1pfW8j+1oAHkkXnePcQg8EJ1kVSRouKgPDv10JFldECoqVvhNeP93DXxtqEClouK89jdVU+T+/xThk9d6CDzHQfV62I/jK3uXJ5eKz9IvWSrVMv8+hfBi5ZGq6a/mYKrQXsauqjKDeD5XMUoS32Z7KszM/Oxr7kCRYDxhj2tQRYV1sw65g1NQUc7RxKydD0vZaiWVdztjKyn2mvKiPlXODJt1swBu66KPrC/2zcsa6a7Sd76Ah4E1X3woEOLl9WOmcU3UzK87NYUZGXEsrozRO9LCv3U5qXFXVsRX42S0pz2Z5CaxhvN/WzvrYw6peZjYuL2NnYm1Kh6W2BUXqGxiJaFjYXVBcwGTIcbk+9KMa9zf1kpvtoiBA4kp+dQXVhNkdTMPrSKVQZLSAe39nMxYuLWFLqT3iOd62vwhg8ycc43jXEsa4hbowhcGEmVywrZfvxHiYmvVt/McbwZmMvm+LoE3XJkmLePJkaL/XR8UkOtQ+wYdHsL3ObdTWFdA2O0TkYTIJksbG3Oex+W1szh2VUHT63L4Vcozb7WgJcUJU/azrGioo8DqsyUlKdd1oDHGgb4K45kkRjYUVFHotLcnnBg2oM9j1vXF0Z97WXLytlaGzSU5/60c4h+obH2bQk9mKcFy0upntojOY+74JGbPa3BpgIman1iblYXZUPwMG21AmT3tfSj0jY+pmNxSW5ZGf4UvKlfrRzkBUV+bOeX1GRx9HOwZQI1HEDVUYLhMd3NZPuE949R6JlLIgIN6wq59WjXUn3q792tJslpbkJVeG2Q3m9jJR682TY3XZxHJbRuhp7LcD7RfU9Vq28WCyjVSmojPa3BKgv85ObObuL1+cT6svyOJZihUcHRsdpDwRZXjG7V6OhIp/hsUlaPIx2dRNVRguAUMjw810tXLuyPKa1imhcb5V8eT2JfY4mQ4bXj3cnFAUI4XWjupIc3mr0Uhn1Wov/sbtJV1cV4BPYnwJuoz3N/ZTlZVIVQ35aaV4W5flZHEghZXSkc5CVc1gWNsvK/RxLsdI6xzrD8swVOLKiInwuFat2OIEqowXA9hM9tPSPxlSHLhauWFZKdoYvqa66/S0BBkYn5sxvicZFdcW8dbLPOaHiZMfJXi5ZXBxXJGNOZhrLy/PY1+K9ZXSwbYALqgtiln91VX7KWEZjEyFOdg9PvbDnYnmZn6ae4ZSq9m63iJhLGdlJ1E09w0mRKdmoMloAPL6rhdzMNG5ZE/9aSySyM9K4cnlZUttjbz0WLi9zxfLEldHFi4toC4x6krTbOzTGsc6hqXDteFhXW+h5rtFkyHCofYCVldEtC5tVlfkcah9gMgXWME50DzEZMjEpo2XleYQMNHanzkv9aOcg6T5hyRwu6vL8LLLSfTSqMlJSkbGJEE/taeXWNZVz+srj5fpV5ZzoHk5apeCtR7tZVuaPq4TRTC5aHFYEXlhHtmVzYQyL/zNZW1NAeyBIl4eRaSe7hwhOhKbWgmJhZVU+wYlQSrwcbddVbMoo7EY92pk6rrqjHUMsLs2ds7CxiLC4JDcl/r3dQJXROc5LhzrpHxnnzo3zi6KbybVW1eBXDrtfO21iMsT2E71cPg+rCMJRVFnpPnZ6sG60vzVs2ayZI5JrNtbU2OHG3rnqDlnFQ1fHoYyWlYVf6idSYP3FVkbLYlivq7fkPtaVOmsvx7oGWVYWXZGGlZEGMCgpyOO7minxZ8ZUfiYelpTmUleSw8tJqM68tyXAYHAi4eAFm8x0H+trCz0JYtjXEqCmMJtif2bc166tDkev7fdQGR1oG0AkHLEVK/ZLPRX67BzpGKS2KCcm70B+dgYV+VlTQQNeY4yhsWd4zsK6NnUluTT1DKdEXprTqDJKMqPjkwRGxx2Za2B0nGf3t/OeDdUx9S2KBxHh2oZyth7tZtzlRNKtR8OVE+YTvGBz8ZJi9rYEkr44vb8lMGXhxEthbgaVBVkc9rC1wcG2AZaW+iPWdJuNEn8mBdnpKaOMYnHR2Sy2XuqpQOdgkNHxUEwpDYtLchkMTtA77Mw7JJVQZZREBoMTXPF3z/Hgy8ccme/X+9oJToQcd9HZXNNQzmBwgl1Nfa7Mb7P1WDcNFXmU588/LP2iuiLGJkJJtTJGxiY52jnImgg1xWJlZWU+hzq8VUaR+hfNhYhQX57nuTIKhQzHuuJTRnUluZzqTQ13l60U64pjU0bAglw38kQZiUiJiGwRkcPWz4ghSCJyrzXmsIjcO+34JSKyR0SOiMi3xYpFnW1eCfNta/xuEbl42lyTIrLL+jzp5nPnZaVzYV0R/7mjyZGyNY/vaqauJIeLFxfNX7gIXLG8lDSfuNphcnwyxI4TPY5YRRCumQaw+1TyotMOtg8QMomtF9k0VORzpMOb7PrR8UlOdA+xqip++ZeV+T1XRq2BUUbHQzGtF9nUFefQ2j/iutUfC03WGlBdSU7Usbb1pMrIOe4DnjPGNADPWftnICIlwP3AZcBm4P5pSuu7wCeBButze5R575g29lPW9TYjxpiN1ud9zj1iZH7/siW0B4I8N88cns6BIK8e6eLOCxOr0B0LhTkZbKwrcnXdaPepfobHJucV0j2dqoJsSv2ZSa1ubNc5m6smWjRWVuYxOh6iqTf5L5kjHYOETHzBCzb1ZX6a+0Y8rYJth2gvKYldGS0qySVkoLXPm4LA07Eto0UxWEaLisMK65QHvydu45UyuhN42Np+GLgrwpjbgC3GmB5jTC+wBbhdRKqBAmPMNhNexXtk2vWzzXsn8IgJsw0osuZJOjesKqeqIJsfvd44r3l++tYpQoZ516KLxjUNZew+1UffsDtNvexK205ZRiJi5e0kz023vyVAQXb61IsiERqs/B4vqkkfttyD8brpAJZaQQwnPczZsV/mtgsrFmyXmBfKfyaNPcNU5GfF1EwyNzOdwpwM2vq9V6JO45UyqjTG2P2t24BI2Zq1QNO0/VPWsVpre+bxueadbS6AbBHZISLbROSuuYQWkU9ZY3d0dibmukpP8/HBS+t4+XBnwguooZDhJ280srm+JC4/eSJc01COMfDqEXfaM2w92s3qqnxKEohCm411tQUcbh9I2rf1fVbwwnwsVLttgBfrRsc7h/AJLI7DsrBZNhVR512YdGPPMGk+oboo9hw1+4tDKgQxNPUOUxeHIq0uzKYlBSw6p3FNGYnIsyKyN8LnzunjLOvGcUd5HPMuMcZsAn4P+CcRWT7HnA8YYzYZYzaVl8ff5sDmQ5vrEODHbyRmHW091s2J7mF+b/PihGWIlQsXFVKQnc5vXMg3Ck5MsuOkc+tFNutqCpmwKgq4zWTIcKAtwJrqxIMXAAqsfjVeWEbHu4dZVJxLZnr8rwPbMvIygbSpd5iaouy4IkqrC7NJ80lKWEZNPSNxWXVVhdm0BVIj+MJJXFNGxpibjTHrInyeANptN5n1M9ICSjNQN21/kXWs2dqeeZw55p1tLowx9s9jwIvARQk+csxUF+Zw4+oKHt1xirGJ+BdQH9l6gqLcjLhacydKelq46+pvDnc5ntvwdlM/o+Mhx9aLbOzmasmohH28a5DR8dC81otsGqzyOsnmeNfgVM5QvORlpVOWl+mphdHYMxzXyxzCv9c1RdlTwQNeMT4ZorV/hLo4XLzVhTkpsdblNF656Z4E7Oi4e4EnIox5BrhVRIqtwIVbgWcsN1xARC63oujumXb9bPM+CdxjRdVdDvQbY1qtubMARKQMuArY7+iTzsI9VyylazDIY2+dij54GofbB3hmXzsfuXxJTD5mJ7imoZzmvhHHKx1vPdqNCFxWH3v/n1hYVJxDQXZ6Uuq92VUTEs0xms7KijyOdAwmtdabMYbjnUMJKyOwEjE9tDCaElBGEF438toyaukbIWTCARWxUl2YTffQWEq2Tp8PXimjrwO3iMhh4GZrHxHZJCIPAhhjeoCvAdutz1etYwCfBR4EjgBHgafnmhd4Cjhmjf+edT3ABcAOEXkbeAH4ujEmKcromoYyNiwq5F9fPBpXmPd3XzpKTkYaH7uq3kXpzuQaq7qD0yHe2451c0FVAUW5zq0Xwekghn1JiKjb3xIgM83nyNrdyspwrbdkRkp1DgQZGpucnzIq9q5e2lBwgq7BsbjWXGzqinM9t4ymwrpjiKSzqS4Mr411BFKny64TeKKMjDHdxpibjDENljuvxzq+wxjziWnjHjLGrLA+3592fIfl8ltujPm8tT4017zGGPM5a/x6Y8wO6/hr1v6F1s//SNa/gYjwhRsbaOwZ5mc7m6NfQPgb4BO7Wrh782JHF/yjUVeSS32Zn984GOI9Oj7Jm429jrvobNbVFvJO24DreST7WwOsrMpzpAKGnSeTzDI1trU7H2W0uCSXlr5RT1q+25ZNQpZRSQ5dg0FPLQy7wnxtUXxuOmDBNdnTCgwecvMFFVy4qJC///VBhoITUcf/45ZDpPmET16bPKvI5pqGMrYe7XaszM7Oxj7GJkLzrkc3G2trChibCLnaiMwYw76WwFRtufmyzOplk8zGbyccUEZ1JTlMhgytHoQb2zlGiSgj+6Xuhdw2doh2RUHs1UfsqMGFFt6tyshDRISvvHct7YEg//bS0TnH7jnVz2M7m/mDq+qn/oiSyTUN5YyMTzrWnuG1o12k+YTNy5xdL7JZa5XmcbMSdnsgSM/QmCPrRQDFuRkU5mQktSX28a4hMtN91MTxzXwmtovMiyCGxgRyjGxsd5cX/a9sWvpHKfVnxrX+a8utlpHiKJcsKeaujTX820tHZ60aEJyY5H/+dDdleZl89oZZI89d5YrlpaT7xLEQ71eOdLFhUSEF2RmOzDeT+jI/ORlprlZisCsvOKWMRIRl5cktr3Osa4glJbmk+RLPkbLXO7xYN2rqGSY/O5wIGi/VlgL20sJo6x+hqjC+Hl65mekUZKerZaQ4z/3vXUupP4sv/ngnvUNnVjoIhQx/9t+7eac1wNd/e4NrL+9o5GWlc/GSYkfWjQKj47zd1MfVK5xtezGdNJ9wQXW+qwVT7bkvmEdNupnUJ7nW2/Gu+UXSQfiberpHOTuNPcPUFecmlHBcVWBbRt691Fv7RxPydJTnZ3najNENVBmlAMX+TL5990Wc6hvh7u9t44iVhT8YnOB//NfbPL6rhT+9bRU3O9RWPFGubShjb0s/3fP8I3j9WA8hA1e5qIwg7Krb3xpwrfjovpYAS0tzyctyrsPu8vI8WvtHGR6LvoY4XyZDhsbuYerjKDAaiXDOTo4nTd8SyTGyyclMoyg3w1M3XVtgdMrtFg8V+dkaTae4w+b6Eh68ZxNtgVFu+ceXuf2fXubKv3uOn+1s5k9uXclnr/fGPTcduzTQK0fmZx29eqSLnIw0LnKp2rjNutoCBoMTrrmP9rcGptamnCKZDeta+kYYmwxNlfSZD170BwqFDE29IzH1AZqN6sIcz9xdI2OT9A2Px+2mg7Bl1KmWkeIW164sZ8sfXceXbmqgujCbW9dW8bPPXsnnb2xwrTJ3PKyrLaQoN2PerrpXjnRxaX0JWenuJu3aisKN5NfA6DiNPcOOrRfZJDO8+3RY9/xzpOpKcpKujDoGgoxNhBLKMbKpLsz2zE3XFhidkiFeyvOz6AgEF1THV+f8C4ojlOdn8eWbV3otRkTSfGKVBurEGJOQgmzuG+FIxyAf3FQXffA8aajMI90n7GsJ8J4NNY7O/Y6DlRems7TUj0hyLKPjVtTe0rLEX+Y2i4pz6R4aYyg4gd9Bt+VczCfHyKaqMNv15pGzYbsHE7GMKvKzGBmfZGhs0lE3sZeoZaTExbUNZbQHghxOMH/nBauP0w2rEy80GytZ6Wk0VOa7Et5tz7nWweAFgOyMNGoKc5IS3n28a4i8rHTK8+bfYddWCMkMYrBzjOKp6zaT6oJsejwqrWO7BxMNYIBwBY2FgiojJS6uaQgrkURLA71woIO6khyWl7vb+sJmXU0B+5r7HXdn7G8NUJaXRUVB/N9qo5Gs8O7j3cPUl/kdcQHbLRla+pIXDGDfaz45UnZ4d3sg+a462z1YlcDvkK2MOjyQ2y1UGSlxUVOUw4qKvITWjUbHJ3n1aBc3rqpI2hrY2poCuofGaHc48sjuYeQGy8r8HOsccn09YD7VumdiK4TmJFaTbukfoSwvvoTRmZxOfPVCGY1QlJtBTmb88lfkh+VeSEEMqoyUuLluZTlbj3UzMDoe13Vbj3YzOh7ihtUVLkl2Nmtr7UoMzgUxhMsMDTjSNiIS9WV+BqwCoG4RnJjkVO/IVD+i+VKel0VGmiTVMmruG52XVQSn12u8CO9u6x9NyCoCddMpCgDvWl/F2ESI5w9EakM1O88f6CAnI83xZnpzcUF1ASLOlgU61D7A+KRhjcPrRTZ2jTo3XXWN3cMYgyNh3QA+n1BVmJ10N13NPEtjeWsZJZZjBFCUk0G6T+hQZaScz1xUV0x1YTa/2N0afbBFKGR4Zl8bVzeUJa0PE4QrR9SX+h21jPa3WsELLlpGgKtBDMcdKJA6k5rCnKQpI2NMWBnN0zLKzUwnPyvdkwTStv7RqTWrePH5JJxrpMpIOZ/x+YQ71lXz0sHOmF1120/00DEQ5D0bql2W7mzW1BQ42vV1f0uA3Mw0lpQ69yKfTk1RDpnpPlctI3tup9x0EG6D0JKkNaP+kXGGxyapKZp/AEl5QRYdA8m1jMYmQnQPjVGZn7j8qowUBXj3hmrGJkM8+057TON/sbuV7AwfN1+Q/JJGa2sKae4boW/YmTWYvc39rK0pmFdx0blI8wlLS3NdbSVxvGuIUn9mQgVGZ6OmKIe2wGhSOtU298XfB2g2KqwE0mRiBx7E0zpiJuV5WeqmU5SL6oqoLcrhp29GbwwYnJjkl3tauWl1ZdISIqezrjbsTnOiaOpkyLhSBmgmbhdMPeZAgdSZ1BSF+xolw8qwLbD5uukAKguyk/5Sty2aivx5KCO1jBQl7Kr70KV1vHKkK+pL8+k9bfQMjXH35sVJku5MnOxtdLxriOGxSdfWi2zqy/I42T3kmpXhRLXumdhN35KxbuREjpFNRX4W7YHRpJbWsfODKubppuseCibFEk0GqoyUhPng5jrSfcL/3XZyznE/3HaS+jI/V7rUYjwaJf5MqguzHalRZwdCrKt11zJaVu5nfNLQ3Ov8i30wOEHnQHDe1bpnUpvEXKOWvhEy032U+jPnPVdlQTbBiRCBUfcrpdvYllj5PCyjEn8mxuCY+9lrVBkpCVORn8271lfz4zca6RmK/AfxxvEedpzs5SOXL8Hn0hpLLKytKXDEMtrXEiAz3ceKCncrSNgh18e6nI+os1uNOxXWbTPVgTQJllFz3wg1hdmO/E55Uc2gYyCICJTlJa5MS60yTrP97Z1rqDJS5sUXb1rB8Pgk//5y5Lbp//zcIcrysjxz0dmsrSnkWOcgQ8H5ffvd29zPBVX5ZKS5+6fjZiuJYy5E0gHkZ2dQkJ2eNDedEy46OO0qS+a6UedAkFJ/Junz+D2yrcJuVUaJIyIlIrJFRA5bP4tnGXevNeawiNw77fglIrJHRI6IyLfFqi0z27wislpEtopIUET+ZMY9bheRg9Zc97n53AuRFRX53LWxlu+/eoKjM/Jint3fzqtHuvn0dcsSKnniJBsXFxEysPtU4q46Ywx7m/tZ43LwAoRdMAXZ6a60kjhuzbnUhdD0mqLk5Bq1OFB9wabSimhLZn26zoFRyuexXgTh3xGAbhcrdSQTryyj+4DnjDENwHPW/hmISAlwP3AZsBm4f5rS+i7wSaDB+tweZd4e4IvA38+4RxrwHeAOYA1wt4iscegZzxv+/F2ryclI44s/3jmVd9TUM8yf/2wPq6vyueeKpd4KCGxcVATAzqbehOc41TtCYHRiKjrPTUSE+vI8VyyjE91D1BbluJJ8XFuU4/qa0fhkiPYB55SRXew2mZZRx0BwXutFcNoy6hlaGBF1XimjO4GHre2HgbsijLkN2GKM6THG9AJbgNtFpBooMMZsM+Hwl0emXR9xXmNMhzFmOzAzQ3MzcMQYc8wYMwb8xJpDiYOK/Gz+6UMbOdg2wHv/5RX+6sl9/Na/vsbYRIh//tBFZKZ77w0u9mdSX+ZnV2NfwnPYfW821BY5IlM0lrkU3u1GWLdNdZH7JYHa+kcxBmodSHiFcJUOf2ZaUi2jjkBwXmHdEP6dBnXTzZdKY4xdS6YNiJQJWQs0Tds/ZR2rtbZnHo913ljuERER+ZSI7BCRHZ2dibVQWKjcsKqCRz6+mdzMdH78RiMNFXk8+odXsKoq32vRpthYV8TOpr6EQ3jfauwlO8PH6urkPFN9mZ/mvhFHe+0YYzjeOehIQ71I1BTl0D8yPu+1ublwMqzbpiKJuUahkKFrcP7KKCPNR2FOxoIJYHAtA1FEngWqIpz6i+k7xhgjIo4HyrsxrzHmAeABgE2bNi2M4H4HuXJ5GU996RqvxZiVixYX8bOdzbT0jyaUub+zsY8NtUWuBy/Y2NbLie4hVlc54xrsGRojMDrhynoRnA7vbu0fYUWFO0q7pd95ZVSen0Vnkqow9A6PMREy81ZGEHbVqWUUBWPMzcaYdRE+TwDtlrsN62ek8s/NwPTe1IusY83W9szjxDhvLPdQFiAb64oA2NkY/7pRcGKS/S0BLlpS5KxQczAVUedgEIPt9nOruWEy+hpNVV+YZ8Xu6VQWZNOepPp0p3OM5u9mLPFn0r1Aehp55aZ7ErCj4+4Fnogw5hngVhEptgIXbgWesdxwARG53Iqiu2fa9bHMO53tQIOI1ItIJvAhaw5lAbK6qoCsdF9C60Z7mwOMTYa4qC5i4KcrTFXvdnDdyI7OW+ZwwquNrYzcXDdq7huhxJ/paISmXZ8uGVUYbGU0n7p0NiX+zAXjpvNKGX0duEVEDgM3W/uIyCYReRDAGNMDfI2wwtgOfNU6BvBZ4EHgCHAUeDrKvFUicgr4Y+AvReSUiBQYYyaAzxNWfO8Ajxpj9rn76IpXZKb72LCokO0n47eMbGvq4sVFDks1O/6sdCoLshwNYjjaNUhGmjhSYDQSlflZ+MRdZRTOMXK23XtlQRYj45MMuLjWZeNEXTqb0rysBaOMkl+1EjDGdAM3RTi+A/jEtP2HgIdmGbcujnnbONO1N/3cU8BTcYivnMNcVl/Kd186ymBwgrw4irbubOyjtihnKgw4WdSX+R3ta3S8c4glpf55JVvORXqaj6qC7Kmq2m7Q3DvieDTgVOJrIEhBtnOVzCNhF5Kdb2g3hNeMeobGCIWMpxVOnCDqb6SIXCEi3xGR3SLSKSKNIvKUiHxORNzP/lMUB7lieSmTIcP2Ez3RB1sYY3irsZeLkmgV2dSXOZtrdKxryPEyQDOpdjHx1ammejOxXWbJKAnUEQiSl5VObub8bYESfyYhA30jsfUVS2XmVEYi8jRhS+UZwoml1YSTQ/8SyAaeEJH3uS2kojjFxYuLyUgTth3rjvmaxp5hWvtHuay+xEXJIrOszE/v8Di9DrhiJiZDnOwecrxA6kxqXGyyFxiZYGhs0nE3YzJLAnUOzD+s26Y0b+EkvkZTzR8xxnTNODYIvGV9viUiZa5IpigukJOZxkV1xWw7Grsy2mqNvcKDquNTEXXdQ1NJjonS3DfC+KRheZm7RV5rirJ5Zu+oK66jZhdyjOC0ZZSMxNdOB6ov2JT6w/N0D46xosKRKT1jTssogiJCRG4SkfeKSMZsYxQllbl8WQl7mvtjbpm+7Vg3ZXlZroVDz4VtxTgR3u12JJ1NbVEOY5MhulwIOT6d8Ors2l1+VjrZGb6kNKvrGBh1bO2xZKok0LkfxBDXKqaIfAu4CriQ6GHTipKSXL68lJAhpnUjYwxbj3Vz+bISrHq8SaWuOJc0nziybmSHiLtVCsjGzv9xI4jBTnh12k0nIpTnJ6eNd4cLbrquha6MRORbIlI07dBiwuHWf2NtK8o5x8WLi8nNTOOFA9FLOh1oG6A9EOSqFd54ozPTfdQV5zijjDoHKczJmPo27Ra1xe4po+a+ETLTfJTlOfMyn05FfrbrltFgcILhsUnH3HTFuZZltAAqd0ezjB4DfiIiX7QqXD8CvABsBb7ntnCK4gbZGWlct7KcX+9vIxSlZfOW/e2IwE0XeOeQry/zO5L4eqxziGXlftctPDcTX1v6Rqkucqap3kwq8rOmwq7dwskcIwh/WSnITl8QAQzR1oxeNcbcTrgFwzOAGGOuN8Zcboz556RIqCgucOvaStoDQXY3z93faMv+djbWFU1FW3lBfVkeJ7qGoirOaBx3sVr3dApzMsjPSncloq6lb8TRMkDTKc/Pct0yskPHnfx9Ks3LWhD16aK56dJF5N2Ea7zdBVwoIk+KyIXJEE5R3OLGVZWk+4Sn9rTOOuZU7zB7mvu5dU2ker/JY1m5n5HxSdrmEek1FJygLTCatCCMmqIcTvW6YRk5n2NkU5GfRWB0wtEq6TNxshSQTbg+3QJXRsDjwEbgOuA7xpivAZ8GviAi6qZTzlkKczO4cXUFj711irGJUMQxP30zXDP3PRuqkynaWSxzoAX5kY5wFYflLkfS2dQWO5/4Oj4Zoj0w6lgfo5nY6zhuWkdTRVIdXPNaKPXpoimjJcaYvwH+GlgPYIxpMcZ8gnCHVEU5Z7n7ssV0DY7x7DvtZ52bDBke3dHE1SvKqCtxp/dPrNjh3TPbusfDwfYBAFY51IoiGjVFzpcEag+MEjLO5xjZJCPxtXMgSGaaj6Jc50oOleUtjDYS0ZTRAyKyFXgJ+IfpJ4wxu9wSSlGSwbUN5dQW5fD9V4+fVa35pUMdNPeN8LuX1s1ydfKoKsgmPyudQ5ZCSYTD7QNkpftYnCTFajfZG3Sw8OhU6wiXlFFyLKNRyvOzHA0iKfFn0js8Nu81Ra+JFsDwL8aYK6zPD5MllKIkgzSf8IfXLWP7iV5ePXK6IoMxhn9+9jC1RTncvtbb9SII58A0VOZxqH0+ltEgDZV5pCWpmGatCxF1zX3DgJuWka2M3Iuoc7L6gk1xbiaTIcPAqPsVx5/a08o/bjnkytzRAhj+0uolNNv5G0XkPc6LpSjJ4Xc31bGoOIe/fHzPVKvsH7/RxNun+vniTSvITPeqy8qZrKrK51D7QML9dg61DbDSpc6rkagtcj7X6LRl5M6aUWleuP2Fq5ZRwHllZCe+dichvPv5Ax08uqPJlbmj1abbA/xCREYJ16LrJFwgtYFwYMOzwN+6IpmiJIHsjDT+/gMX8nvf28bvfPc11tYU8viuZq5pKOMDl3jvorNpqMjnx8NNdA4G4w4L7h8epy0wysqqJCqjYjcso3BTPSeqXUcizSeU+N2twtA5GGTTUmcbNNqJr73D7q8b9Q6NuZY0Hc1N94Qx5irCEXT7gDQgAPwQ2GyM+SNjTPQ0dkVJYS5fVsq/f2QT45Mhtuxv4/0XL+LfPnxJSvWHWVkZViSHE3DVHeqwghcqk6eMKvKzSfcJzQ6Gd7vRVG8mFS6WBBqbCNEzNOZ4ztr0Yqlu0zPsnjKK6SuGMeYwcNgVCRQlBbhlTSW3rKn0WoxZWVkVzg861D4Qd2mig20D1hzJU0ZpPqGqMNtRy6ilb4Slpe6GpruZ+GoXjnV8zcgfjsxLlmVUV+xOEExqOMQVRZmT8rwsinIzEoqoO9Q+QF5WOjWFya0iUVOU49iakTGG5l73El5t3CwJ5HQpIJspyygJ4d3dXrnpFEVJDUSElRX5CUXUHWwbYGVlXtKrji9ysMleYNSdpnozKc/PomvQnTBpWxk5bRnlZKaRneFzpAHjXIxPhhgYnZhao3KaWNqOp4nIH7lyd0VRYmZlVV7cEXXGGA61D7AqiS46m5qiHNoCo0xMRq5wEQ8tLjXVm0lFfhaTIUOPCy6vzkHnSwHZlPrdr09nuwFL/M4l7E4nqjIyxkwCd7tyd0VRYmZVZT4DoxO09MdubXQOBOkdHqchiWHdNjVFOUyGDO0OrMG41VRvJuVWcIEb60YdgfCctlvNSYr9Ga5bRr1D4WaUJS7ID7G76V4Vkf8tIteIyMX2J9GbikiJiGwRkcPWz4ixjiJyrzXmsIjcO+34JSKyR0SOiMi3xfI/zDaviKwWka0iEhSRP5lxjxPWXLtEZEeiz6QobrO2thCAvVEqjU9nX0sgfG1NcsoATcfJ8G57DntOt7CtFjci6joHRynOzXAld63En+V6fTp7/mKvLCOLjcBa4KvAt6zP38/jvvcBzxljGoDnrP0zEJES4H7gMmAzcP80pfVd4JOE850agNujzNsDfHEOmW8wxmw0xmyaxzMpiqusqS4gzSdxKSN7rK3Ikold0NSJ8O5TdlM9l76V21S4WBKocyD+HLFYKcnNcMW1OJ3TbjpvQ7tvcPi+dwLXW9sPAy8CfzZjzG3AFmNMD4CIbAFuF5EXgQJjzDbr+COE21s8Pdu8xpgOoMNqh6Eo5yTZGWmsKM9jTzzKqKWfZWV+8rLcSRSdixoHqzC42VRvOnZwgRsRdR0ulAKyKfFnTbnR3MJekyrxKoABQEQKReQfRGSH9fmWiMznq1alMcZuJNMGRErwqAWm1504ZR2rtbZnHo913pkY4Nci8qaIfCpG+RXFE9bVFrK3uT/mIIa9zQFPrCKA3Mx0inMzHFJG7jXVm05uZjp5WelT6ztO4kZdOpsSfwaDwQmCE+71YrLXpIq8VEbAQ8AA8LvWJwB8f64LRORZEdkb4XPn9HEm/FfleBxlHPNebYy5GLgD+JyIXDvbQBH5lK2QOzu18ISSfNbXFtA1OEZ7DC/L3qExmvtGWOfBepGNU32N3GyqN5Py/KypyDenMMa4rIzC87ppHfUMjZGfle5avcZYbfflxpjfmbb/1yKya64LjDE3z3ZORNpFpNoY0yoi1YQ7yc6kmdMuN4BFhN1uzdb29OPN1nYs886Us9n62SEiPyO8PvXyLGMfAB4A2LRp07ldr105J1m/KGzl7GnupypKEuvbp/rC13hkGQHUFOZwojvxpoDgflO9mZTnZ9HpsGU0EJwgOBFyPOHVxg637h4KRv29SJTe4TFK8tyxiiB2y2hERK62d0TkKmA+X3eeBOzouHuBJyKMeQa4VUSKrcCFW4FnLDdcQEQut6Lo7pl2fSzzTiEifhHJt7ete+xN/LEUxV3WVBeS5hPebuqLOnbHiV7SfMLGxUWuyzUbtcU5NPeOJFxtHNxvqjcTNywj2+3nlmU0VSzVZcvIrYRXiN0y+jTwyLR1ol5Ov/QT4evAoyLyceAkYdcfIrIJ+LQx5hPGmB4R+Rqw3brmq3YwA/BZ4AdADuHAhaejzFsF7AAKgJCIfBlYA5QBP7Miw9OBHxljfjWP51IUV8nJTGNtTQHbT/REHfvGiR7W1RS4VuU6FmqLchgamyQwMkFhgt1N3W6qN5OK/CxeDDgbwNDpQrvx6SSjjUTvsPNFXqcT9bdURNKAjxhjLhSRAgBjTGA+NzXGdAM3RTi+A/jEtP2HCK9XRRq3Lo552zjTtWcTAC6MR3ZF8ZrNS0t4ZNtJghOTZKWnRRwTnJjk7aY+Pnz5kiRLdya2AjnVN0xhbmLuwlO94aZ6bucY2ZTnZzE0NslQcAK/Q1GInS4VSbU5bRm5F97dMzjGqkr31h9jrcBwtbUdmK8iUhRlfmyuL2FsIsTuU7OHeO9t7ic4EeJSh3vnxMvpjq+JWxpNPSOI4HpdOpsKF6owdFiWlluWRVFuJiK4mvgabh/hTsIrxO6m2ykiTwL/BUytRhpjHnNFKkVRZuXSpSUAvH6se2p7Jq8c7kaEWc8ni6lcI8u6SYTGnmEq87PJzohsBTrNVOLrYJClZc60rOgcDJKZ5qMgx73GgEU57iW+joxNMjoeotilhFeIPYAhG+gGbgTea3203biieECxP5MLqgt4+XDXrGOeP9jBhkVFlLq0RhErZXmZ5GSk0diTeLxTU+8wi0vc6aETianEVwcj6uywbjcrp5f4M12zjGwl51bCK8S+ZtRtjPmTaGMVRUkON19QwXdeOEJPhP4yXYNBdp/q48s3rfRIutOICItLcmnsSTy8u6lnmCuWlzoo1dycLgnkXBCDmzlGNm4qI3styq1SQBD7mtFVrkmgKErc3LKmkpCB5w+cnUr30sFOjIEbV1d4INnZLC7N5WR3Ym664MQkbYHRpFpGxbmZpPvE0WKp57oy6kkFZWSxS0SeFJGPiMhv2x/XpFIUZU7W1xZSVZDNM/vazjr35Nst1BRme1KpOxJLSnJp7BlOqGFdOEeJpCojn08oy8s6R5WRO3lGpyt2e6+MdM1IUVIIEeHOjTU8f6CDtmn9jdr6R/nN4U5+55JFrhcVjZUlpbkEJ0IJvdwbe8IWVV0SlRFYia8OKaPxyRA9w2Ou5RjZlPgz6R12p0ttj8tFUiFGZWSM+ViEzx+4JpWiKFH58OVLCBnDw1tPTB17eOsJDPA7F0dKq/OGJaXhiLSTCZQFarLaTyTTMoLwupFTllH34BjGuNPhdTrFuZlMhgwDoxOOz907PIZPoCDHvdDuWKt2rxSR50Rkr7W/QUT+0jWpFEWJSl1JLu/dUMNDrxznndYAB9sG+P6rx3nPhhrHQpKdYElpWJGc7Il/3aipZ5jMdJ/rVsVMnLSM3K6+YONmFYaeoTGKcjNJc9HajtVN9z3gz4FxAGPMbuBDbgmlKEps/MW7L6AoN4O7vvMqd37nFfKy0vlf777Aa7HOoKYohzSf0JhAEENTzzB1xTlJdzlW5GfRPRRkYjI077k6B8NuVLfXjKaqMLiQa9Q7fHbUptPEmoGVa4x5Y0aMvPO2oKIocVFZkM1PP3Ml//bSUSZD8JnrllNRkJzq1rGSkeajtignIcvoeNfQlJsvmZQXZGNM2CKY77+nbRm5/f9SarWRcCOIoWdozNX1IohdGXWJyHKs/kAi8n6gde5LFEVJBouKc/l/71rvtRhzsqQ0l8Y414xCIcOJ7iGuXlHmklSzY7vUOgaC81YidvJsmYvtFwCKrVI9PS656epddv3G6qb7HPDvwGoRaQa+TLiSt6IoSlQWl+RyIk43XUv/CKPjIZaV57kk1ezYwQZOrBt1DgYpzMmYtaitU7hrGY2nhpvOGHMMuNnq+eMzxgy4KpWiKAuKJaW59I+M0z88HnMriWOdYUtqWbkHbropy2j+VRg6Au7nGEG4vUh2hs9xy8gYQ++wu72MIHbLCABjzJAqIkVR4sVe9zkeh6vueJeljDyIDCzPd84yah8YpSpJ63il/izHLaPA6ASTIeO6ZeROM3NFUZRpNFSEXW1HOgZjvuZY5yB5WelJsSpmkp2RRkF2uiO5Ru39o67nGNkU+zMct4ySUZcOVBkpipIEFpfkkpnm43BH7I6VY11DLCv3u1rpei4qCrLnXbk7FDJ0DASTZhmV+LPoGXbWMrIrdrtZCghiT3r9gIjkW9t/KSKPicjFrkqmKMqCIT3Nx9KyXI7GZRkNeeKis6nIz5rq0Joo3UNjTIQMlclSRrnOW0Y9g+6XAoLYLaP/ZYwZEJGrgZuB/wC+655YiqIsNBoq8jkcozIaHZ+kuW+E+rLkR9LZlOdnzTuAod3q8Jo0ZeTPotfhNaOpXkapYBkBk9bPdwMPGGN+CbgrmaIoC4rlFXk09QwzOj4ZdezRzkHrGo8to4EgxiReePS0MkrOmlGJP4PB4ATBiej/xrHSm4SK3RC7MmoWkX8HPgg8JSJZcVyrKIpCQ0UeIXM6ZHsuDraF15ZWV+W7LdaslOdnMToeYiCYeLGZdmvNqaoweZYR4Kh11DM8RmaaD3+mu3lSsSqU3wWeAW4zxvQBJcCfuiWUoigLjxV2RF1ndFfdwbYBMtN9LPWgFJBNRX5YgcwnvLstMIoIlCWp0GuJVYXByWKpvVY3YbcDSWJVRtXAL40xh0XkeuADwBuJ3lRESkRki4gctn4WzzLuXmvMYRG5d9rxS0Rkj4gcEZFvi/WvNNu8IvL7IrLbuuY1Eblw2ly3i8hBa677En0mRVHmpr7Mj0/gSHv0iLp32gZYUZ5Hepp3Dhg7pHw+EXUdgVFK/VlkJOk5XLGMhsZdd9FB7Mrop8CkiKwAHgDqgB/N4773Ac8ZYxqA56z9MxCREuB+4DJgM3D/NKX1XeCTQIP1uT3KvMeB64wx64GvWc+AiKQB3wHuANYAd4vImnk8l6Ios5CdkUZ9mZ/9rdGV0cG2gKcuOgivGQHziqhrC4xSVZi8PCk3LKOeoeDUvG4SqzIKGWMmgN8G/sUY86eEraVEuRN42Np+GLgrwpjbgC3GmB5jTC+wBbhdRKqBAmPMNhNeWXxk2vUR5zXGvGbNAbANsDuPbQaOGGOOGWPGgJ9YcyiK4gLrawvZ09w355i+4THaA0FWeayMTltGiUfUtQeCVOYnr4r6VBuJIefaSPQOj7teCghiV0bjInI3cA/wC+vYfFRlpTHGrvrdBlRGGFMLNE3bP2Udq7W2Zx6Pdd6PA09HuUdERORTIrJDRHZ0dnbONkxRlFlYv6iI9kBwzpDpdyzLyWtlVJiTQWa6b16WUUdglMokBS8AFOVmInK6TbgT9Ay538sIYldGHwOuAP7GGHNcROqB/zPXBSLyrIjsjfA5w/KwrBvHm7ZHmldEbiCsjP4swTkfMMZsMsZsKi8vd0BKRTm/WF9bCMDe5v5Zx+w+1XfGWK8QEcrzsuhMcM0oODFJ99BYUi2jNJ9QlJMxlRs0XyYmQ/SPJMcyirVq937gi9P2jwPfiHLNzbOdE5F2Eak2xrRabreOCMOageun7S8CXrSOL5pxvNnannVeEdkAPAjcYYzpnnaPulnmUhTFYdbWFCACu0/1c+PqSI4L2NXUx+KSXEqT3Go8EuHE18SUkR2Fl8w1IwgnpzoVwNA3Ep6n1OVeTBB7OaAGEflvEdkvIsfszzzu+yRgR8fdCzwRYcwzwK0iUmwFLtwKPGO54QIicrkVRXfPtOsjzisii4HHgI8YYw5Nu8d2oEFE6kUkk3Ar9Sfn8VyKosyBPyud5eV57D41u2W0s7GPjXVFyRNqDuzE10SwE16T3Xm3xJ/pWADDVMJrCq0ZfZ9wBNsEcAPhoIEfzuO+XwduEZHDhMsLfR1ARDaJyIMAxpgewpFv263PV61jAJ8lbOUcAY5yeg0o4rzAV4BS4F9FZJeI7LDuMQF8nrDiewd41Bizbx7PpShKFDYtKWbHiR4mQ2d751v7R2gLjKaMMppPSaCphFcPlJFTllF3kip2Q+xtx3OMMc+JiBhjTgJ/JSJvEn7Jx43lJrspwvEdwCem7T8EPDTLuHVxzPuJ6fPOOPcU8FQc4iuKMg8uX1bKT7Y38U5rgHUz1oVeOdw1NSYVqMjPpnd4nLGJEJnp8eUKtfUnty6dTYk/kzdP9jkyVypaRkER8QGHReTzIvJbgHcVDBVFOWe5bFkJANuOdZ917pUjXZTlZXmeY2Qz1X48gYi69oFRMtN8FMfY2dYpSvyZ9A6Pzaumnk2yiqRC7MroS0Au4SCGS4CPcHptRlEUJWaqC3NYVu7nxYNnpkdMhgyvHuni6hWl+Hze9DCaSbUVlt3aNxL3tXZTvWT3YyrOzWQyZAiMJF5Tz+Z0kdQUSXo1xmw3xgwaY04ZYz5mjPltY8w2t4VTFGVhcvvaKrYe66ZvWgjy68e76Roc46YLIkfZeUFtUQ4AzQkoo+a+EWqs65OJHfnmRHh399AY/sw0stLdLZIKsUfTbRKRn4nIW1aNt90isttt4RRFWZjcvq6KyZDhqT1tU8ee2NmCPzONm1NIGVVbyqS1P/4ghpa+URZ5oIzs9R0nmuz1DI0lLcQ+1gCG/0u4SvceIOSeOIqinA+sry1kbU0BD75yjA9eWkdgZJyf727hXeuryXG5VUE85GWlU5iTQUucltHEZIi2wKg3lpFVLLXHgYi67sGxpOQYQezKqNMYo/k3iqI4gojw2etX8LkfvcW/PH+Yxu5hRsYn+dS1y7wW7SxqinLiVkZtgVEmQ4baYg8sI2t9xwnLqGswyKLi3HnPEwuxKqP7rfyf54CpJzTGPOaKVIqiLHjetb6K29dW8U/PHgbgCzeuYGVlakTRTae2KJvmvvjcdM29I9a157Zl1DM0lrScr1iV0ceA1YSLo9puOkO4qoGiKErciAjf+f2LeXpvK/6sdK5fmZr1HmuKcth+ojf6wGm09FvKyAPLKCczjewM37wto1DIWGtGqeWmu9QYs8pVSRRFOe9I8wnv2VDjtRhzUlOUQ//IOIPBCfKyYntl2pZRTWHylRGEraP5WkaB0XEmQmaqYZ/bxJpn9Jo2nVMU5XzEDkKIJ9eouW+EUn+mZ8EYxf6MeVtGXYPh0PCyFLOMLgd2ichxwmtGQrhLwwbXJFMURUkBaovCia/NfSM0xLim1dw36omLzqbEn0XP8Pwso26r6kRpkiyjWJXR7dGHKIqiLDxsy6gljiCG5t5hGiq8C8Yoyc3geNfgvOawG/Sl1JqRVRxVURTlvKMiP5t0n9DcNxzTeGMMzX0jXL+qwmXJZqfEnzXvyt1dSVZG8ZWhVRRFOc9I8wm1xTk09sS2ZtQ9NMboeMiTsG6bEn8Gg8EJghOTCc9hu+lKklCxG1QZKYqiRGVxSS6N3UMxjT3ZHbaglpYlJ1k0EnYE3Hyso+7BMYpyM0hPS46aUGWkKIoShSWluZzojs1Nd6IrrLSWlvrdFGlOSqaqMCReLLV7KEhpElpH2KgyUhRFicLSUj/9I+NnVBmfjZPdQ/iEpJXRiUTJVBWGxJVR12DyiqSCKiNFUZSoLC4JK5aTMVhHJ7qHqS3OibszrJNMWUbzaCPRMzSWtBwjUGWkKIoSlSWWy+1kT3RldLJ7yFMXHUyzjBLoUGvTPRhMWo4RqDJSFEWJypRl1BU9iOFE9zBLSr1z0QEU5mQgQsKJrxOTIXqHx5MW1g2qjBRFUaKSk5lGZUFWVMuob3iM/pFxzy2jNJ9QlJN4SSDbvbfg14xEpEREtojIYetn8Szj7rXGHBaRe6cdv0RE9ojIERH5tlhN5mebV0R+3+pOu0dEXhORC6fNdcI6vktEdrj97IqinJssKfFzMkp497EUiKSzKfFnJhza3W3VpTsfounuA54zxjQQ7pF038wBIlIC3A9cBmwm3FPJVlrfBT4JNFgfu1zRbPMeB64zxqwHvgY8MON2NxhjNhpjNjn0fIqiLDBiCe8+3D4AkBJ9mUr8mXQnahkNnT/K6E7gYWv7YeCuCGNuA7YYY3qMMb3AFuB2EakGCowx24wxBnhk2vUR5zXGvGbNAbANWOTo0yiKsuBZXpFH50CQ/jnWYQ62DZKTkcYiD4uk2szHMuqyi6QudDcdUGmMabW224DKCGNqgaZp+6esY7XW9szjsc77ceDpafsG+LWIvCkin5pLaBH5lIjsEJEdnZ2dcw1VFGWBsaoqbO0caAvMOuZwxwANlXn4fJIssWYlbBklFtrdneT2ERB71e64EZFngaoIp/5i+o4xxoiIcfr+keYVkRsIK6Orpx2+2hjTLCIVwBYROWCMeXmWOR/AcvFt2rTJcZkVRUldLqgqAOBA2wCXLSuNOOZg2wDXNKRGx9oSfya9w2MYY7CW1WOmeyhImk8oyM5wSbqzcU0ZGWNunu2ciLSLSLUxptVyu3VEGNYMXD9tfxHwonV80Yzjzdb2rPOKyAbgQeAOY0z3NDmbrZ8dIvIzwutTEZWRoijnL5UFWRTmZHCgbSDi+b7hMToGgqyqykuyZJEpzs1kMmQIjExQmBufUukcCFKWl5lUC88rN92TgB0ddy/wRIQxzwC3ikixFbhwK/CM5YYLiMjlVhTdPdOujziviCwGHgM+Yow5ZN9ARPwikm9vW/fY69xjKoqyUBARVlXlz+qmO9Qe7h8UawM+t7FzhBKpwtA5EKQ8P3nrReCdMvo6cIuIHAZutvYRkU0i8iCAMaaHcOTbduvzVesYwGcJWzlHgKOcXgOKOC/wFaAU+NcZIdyVwCsi8jbwBvBLY8yvXHpmRVHOcS6oyudQ2wCh0Nle+r3N/QCsqS5ItlgRKbZaPySSa9QxEKQiP9tpkebENTfdXFhuspsiHN8BfGLa/kPAQ7OMWxfHvJ+YPu+048eAC2ceVxRFicQF1QUMjU1yvHuI5eVnuuN2NvVRU5hNZUFyX+KzUTpVLDX+iLqOgSDrawudFmlOtAKDoihKjGxaGk513HGi56xzb53s5aLFEfP3PaF4qo1EfJbRZMjQPXj+uOkURVHOOZaX51Gcm8H2E71nHO8IjNLcN8JFi4u8ESwCtmUUb3h391CQkIEKVUaKoiipiYiwaWkJbxw/0zJ6+XAXAJfVRw759oKczDTystLpHIjPMuoIhMerZaQoipLCXNtQRmPPMEc6Tod4P7u/ncqCLNbVpkbwgk1FQdaUcomVzkFbGSV37UuVkaIoShzcujacy/+rvW0AjIxN8pvDndx0QWXcyaVuU5mfTcfAaFzXdFrKS910iqIoKUxlQTYXLy7i8V0thEKGx3c1MzQ2yV0ba6NfnGQqCrJoT9gyUmWkKIqS0nzkiiUc6RjkB6+d4F+eO8z62kIuXZo6kXQ2lQXZtAdGCdeUjo2OwCgF2elkZ6S5KNnZqDJSFEWJk/dsqOGixUV89Rf76Roa46/vXJtyLjoIu9qCEyECIxMxX9PhQfUF8CjpVVEU5VwmI83H9z96KY+91cympcVsWFTktUgRsRNw2wdGY65P1+lB9QVQy0hRFCUhinIz+YOr61NWEcE0ZRSIPYihYyBIRUHyLSNVRoqiKAsUOyIu1vBuYwwdA6OUJ7Gpno0qI0VRlAWKbeG0xxje3T8yzuh4iKpCddMpiqIoDpGbmU5+dnrMllFLX1hp1RQlv226KiNFUZQFjB3eHQut/SMAahkpiqIozlKRnxWHMrIso0K1jBRFURQHCVtGsbnpWvtHSPeJJ3lGqowURVEWMNWFYTfdZITutDNp7RulsiCbNF/yE3hVGSmKoixgaotzmAiZmFx1rf2jnqwXgSojRVGUBU2tFRnX3DcSdWxr/wjVqowURVEUp1lUHFZGp3qH5xxnjKG1f9STsG5QZaQoirKgqS3KBaC5d27LqGdojOBEiKoCtYwURVEUh8nJTKPUnxnVTdfYE7ac6kpykyHWWXimjESkRES2iMhh62fEZiAicq815rCI3Dvt+CUiskdEjojIt8Wq3z7bvCJyp4jsFpFdIrJDRK6Odg9FUZSFQG1xDqeiWEa2Mlpaep4pI+A+4DljTAPwnLV/BiJSAtwPXAZsBu6fprS+C3wSaLA+t0eZ9zngQmPMRuAPgAdjuIeiKMo5T21RTlTL6GT3eWoZAXcCD1vbDwN3RRhzG7DFGNNjjOkFtgC3i0g1UGCM2WbCLQwfmXZ9xHmNMYPmdLtDP2BvR7yHI0+oKIqSAtQW5dDcOzJnx9eT3cNUFWQnvcOrjZfKqNIY02pttwGVEcbUAk3T9k9Zx2qt7ZnH55xXRH5LRA4AvyRsHc11j7MQkU9ZLr4dnZ2dUR5PURQlNVhcmktwIjRnJYaT3UMs9shFBy4rIxF5VkT2RvjcOX2cZbHE3qQ9RmbOa4z5mTFmNWFr6WsJzPeAMWaTMWZTeXm5c4IqiqK4yLKyPACOdQ7OOuZkzzBLPHLRgcvKyBhzszFmXYTPE0C75W7D+tkRYYpmoG7a/iLrWLO1PfM4scxrjHkZWCYiZXPcQ1EUZUGwrNwPwNGuoYjnh8cm6BwIsmShWkZReBKwI9fuBZ6IMOYZ4FYRKbaCCm4FnrHccAERudyKortn2vUR5xWRFdMi7i4GsoDu2e7h7KMqiqJ4R1VBNjkZabNaRnYk3ZJSfzLFOoN0z+4MXwceFZGPAyeB3wUQkU3Ap40xnzDG9IjI14Dt1jVfNcb0WNufBX4A5ABPW59Z5wV+B7hHRMaBEeCDlhtvrnsoiqKc8/h8wrJyP8c6I1tGRzvCx+vLzkNlZIzpBm6KcHwH8Ilp+w8BD80ybl0c834D+MYsskS8h6IoykJhWXkeu5p6I5472D6AT2BFRV6SpTqNVmBQFEU5D1hW5udU7wij45NnnTvUNsDSUr9nYd2gykhRFOW8YHlFHsYQ0VV3qGOAhkrvrCJQZaQoinJesLamAIC9Lf1nHB8dn+RE1xCrKvO9EGsKVUaKoijnAfWlfvKy0tnbfKYy2t8aIGRgjaWsvEKVkaIoynmAzyesrSlgzwxltLOxD4CLFntbklOVkaIoynnC+tpC9rcEGJsITR3b2dhLbVEOlR71MbJRZaQoinKecNmyUoITIXacCKdSGmN462QvFy0u8lYwVBkpiqKcN1y5vJTMNB8vHgoXej7SMUhL/yhXLC/1WDJVRoqiKOcN/qx0Lq0v5tl32jHG8My+NgBuWh2paUJyUWWkKIpyHnHXxlqOdQ7xzL42fvR6I1cuL6Wq0Nv1IlBlpCiKcl7xvo011Jf5+fQP36I1MMrnb1jhtUiAt4VSFUVRlCSTlZ7GQx+9lO+8cISrV5Rx5Yoyr0UCVBkpiqKcd9SX+fn7D1zotRhnoG46RVEUxXNUGSmKoiieo8pIURRF8RxVRoqiKIrnqDJSFEVRPEeVkaIoiuI5qowURVEUz1FlpCiKoniOGGO8luGcREQ6gZMJXl4GdDkoTqqx0J8P9BkXAgv9+SD1nnGJMaY80glVRh4gIjuMMZu8lsMtFvrzgT7jQmChPx+cW8+objpFURTFc1QZKYqiKJ6jysgbHvBaAJdZ6M8H+owLgYX+fHAOPaOuGSmKoiieo5aRoiiK4jmqjBRFURTPUWWUZEQkTUR2isgvvJbFDUTkhIjsEZFdIrLDa3ncQESKROS/ReSAiLwjIld4LZNTiMgq6//O/gRE5Mtey+U0IvJHIrJPRPaKyI9FJNtrmZxGRL5kPd++c+H/UDu9Jp8vAe8ABV4L4iI3GGNSKdHOaf4Z+JUx5v0ikgnkei2QUxhjDgIbIfzFCWgGfualTE4jIrXAF4E1xpgREXkU+BDwA08FcxARWQd8EtgMjAG/EpFfGGOOeCvZ7KhllEREZBHwbuBBr2VREkNECoFrgf8AMMaMGWP6PBXKPW4CjhpjEq00ksqkAzkikk74y0SLx/I4zQXA68aYYWPMBPAS8NseyzQnqoySyz8B/xMIeSyHmxjg1yLypoh8ymthXKAe6AS+b7lbHxQRv9dCucSHgB97LYTTGGOagb8HGoFWoN8Y82tvpXKcvcA1IlIqIrnAu4A6j2WaE1VGSUJE3gN0GGPe9FoWl7naGHMxcAfwORG51muBHCYduBj4rjHmImAIuM9bkZzHcj++D/gvr2VxGhEpBu4k/MWiBvCLyIe9lcpZjDHvAN8Afg38CtgFTHopUzRUGSWPq4D3icgJ4CfAjSLyQ29Fch7rWyfGmA7Caw2bvZXIcU4Bp4wxr1v7/01YOS007gDeMsa0ey2IC9wMHDfGdBpjxoHHgCs9lslxjDH/YYy5xBhzLdALHPJaprlQZZQkjDF/boxZZIxZStj98bwxZkF9GxMRv4jk29vArYTdBQsGY0wb0CQiq6xDNwH7PRTJLe5mAbroLBqBy0UkV0SE8P/hOx7L5DgiUmH9XEx4vehH3ko0NxpNpzhJJfCz8N836cCPjDG/8lYkV/gC8H8tV9Yx4GMey+Mo1heJW4A/9FoWNzDGvC4i/w28BUwAOzmHyubEwU9FpBQYBz6X6oE2Wg5IURRF8Rx10ymKoiieo8pIURRF8RxVRoqiKIrnqDJSFEVRPEeVkaIoiuI5qowUJQWwKoF/dtp+jRV+7Ma97hKRr8xxfr2I/MCNeyvKbGhot6KkACKyFPiFMWZdEu71GvC+uSqri8izwB8YYxrdlkdRQC0jRUkVvg4st3oIfVNElorIXgAR+aiIPC4iW6x+UZ8XkT+2CrVuE5ESa9xyEfmVVaT2NyKyeuZNRGQlELQVkYh8wOp587aIvDxt6M8JVwpRlKSgykhRUoP7CLdr2GiM+dMI59cRLulyKfA3wLBVqHUrcI815gHgC8aYS4A/Af41wjxXEa48YPMV4DZjzIWEC6Pa7ACumcfzKEpcaDkgRTk3eMEYMwAMiEg/YcsFYA+wQUTyCBf7/C+rHBNAVoR5qgm3wLB5FfiB1WDusWnHOwhXtFaUpKDKSFHODYLTtkPT9kOE/459QJ8xZmOUeUaAQnvHGPNpEbmMcNPHN0XkEmNMN5BtjVWUpKBuOkVJDQaA/EQvNsYEgOMi8gEACXNhhKHvACvsHRFZbox53RjzFcIWk92AbSULrOK6ktqoMlKUFMCyRl61ggm+meA0vw98XETeBvYRbiA3k5eBi+S0L++bIrLHCpZ4DXjbOn4D8MsE5VCUuNHQbkU5zxCRfwZ+box5dpbzWcBLhLv2TiRVOOW8RS0jRTn/+Fsgd47zi4H7VBEpyUQtI0VRFMVz1DJSFEVRPEeVkaIoiuI5qowURVEUz1FlpCiKoniOKiNFURTFc/5/O+tuUHvuACQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot wave elevation of a standing wave\n",
    "mas = np.loadtxt(\"mass.dat\", comments=\"#\", delimiter=\" \", unpack=False)\n",
    "plt.figure(1)\n",
    "plt.plot(mas[:,0],mas[:,1])\n",
    "plt.xlabel('time (s)')\n",
    "plt.ylabel('mass error (%)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f831a4dc-e99b-4a14-83fd-7f2022017828",
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (131033156.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"/tmp/ipykernel_3928785/131033156.py\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    ffmpeg -i C%4d.png -v -8 C.gif\u001b[0m\n\u001b[0m              ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "ffmpeg -i C%4d.png -v -8 C.gif"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python(Proteus Dev)",
   "language": "python",
   "name": "python-proteus-dev"
  },
  "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}