{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Test of Structural Mechanics\n", "In this example a beam, consisting of two elements, three nodes, is created.\n", "The left most node is fixed and a force is applied at the right most node.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import akantu as pyaka\n", " \n", "import copy\n", "import numpy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating the Mesh" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Create a mesh for the two dimensional case\n", "el_type = pyaka._bernoulli_beam_3\n", "beam = pyaka.Mesh(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now create the connectivity array for the beam." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "beam.addConnectivityType(el_type)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need a `MeshAccessor` in order to change the size of the mesh entities." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "beamAcc = pyaka.MeshAccessor(beam)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we create the array to store the nodes and the connectivities and give them their size. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "nb_elem = 30\n", "L = 2\n", "beamAcc.resizeConnectivity(nb_elem, el_type)\n", "beamAcc.resizeNodes(nb_elem + 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Setting the Nodes" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "Nodes = beam.getNodes()\n", "length = L / nb_elem\n", "for n in range(nb_elem + 1):\n", " Nodes[n, :] = [n * length, 0., 0.];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Setting the Connections" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "Conn = beam.getConnectivity(el_type)\n", "\n", "for e in range(nb_elem):\n", " Conn[e, :] = [e, e + 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Ready\n", "We have to make the mesh ready." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "beamAcc.makeReady()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating the Model" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "model = pyaka.StructuralMechanicsModel(beam)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "normal = beam.getDataReal(\"extra_normal\", el_type)\n", "\n", "for e in range(nb_elem):\n", " normal[e, :] = [0, 0, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Setting up the Modell" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Creating and Inserting the Materials" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat1 = pyaka.StructuralMaterial()\n", "mat1.E = 1e9\n", "mat1.rho = 10.\n", "mat1.I = 1.\n", "mat1.Iz = 1.\n", "mat1.Iy = 1.\n", "mat1.A = 1.\n", "mat1.GJ = 1.\n", "model.addMaterial(mat1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Initializing the Model" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "model.initFull(pyaka.AnalysisMethod._implicit_dynamic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Assigning the Materials" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "materials = model.getElementMaterialMap(pyaka.ElementType._bernoulli_beam_3)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "materials[:, :] = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Setting Boundaries" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "16\n" ] } ], "source": [ "# Neumann\n", "F = 1e4\n", "no_print = int(nb_elem / 2 + 1)\n", "\n", "print(no_print)\n", "\n", "# Apply a force of `10` at the last (right most) node.\n", "forces = model.getExternalForce()\n", "forces[:, :] = 0\n", "forces[no_print, 1] = F" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Dirichlets\n", "# Block all dofs of the first node, since it is fixed.\n", "# All other nodes have no restrictions\n", "boundary = model.getBlockedDOFs()\n", "boundary[:, :] = False\n", "\n", "boundary[0, 0] = True\n", "boundary[0, 1] = True\n", "boundary[0, 2] = True\n", "\n", "boundary[nb_elem, 1] = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the System" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Set up the system\n", "deltaT = 1e-6\n", "model.setTimeStep(deltaT)\n", "solver = model.getNonLinearSolver()\n", "solver.set(\"max_iterations\", 100)\n", "solver.set(\"threshold\", 1e-8)\n", "solver.set(\"convergence_type\", pyaka.SolveConvergenceCriteria.solution)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "model.assembleMatrix(\"M\")\n", "M = pyaka.AkantuSparseMatrix(model.getDOFManager().getMatrix(\"M\"))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "model.assembleMatrix(\"K\")\n", "K = pyaka.AkantuSparseMatrix(model.getDOFManager().getMatrix(\"K\"))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def Mmat2D(L=1, A=1, rho=1):\n", " return rho*A*L/420 * np.matrix([\n", " [140, 0, 0, 70, 0, 0],\n", " [ 0, 156, 22*L, 0, 54, -13*L],\n", " [ 0, 22*L, 4*L**2, 0, 13*L, -3*L**2],\n", " [ 70, 0, 0, 140, 0, 0],\n", " [ 0, 54, 13*L, 0, 156, -22*L],\n", " [ 0,-13*L, -3*L**2, 0, -22*L, 4*L**2]])" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def Kmat2D(L=1, A=1, E=1, I=1):\n", " return np.matrix([\n", " [ E*A/L, 0, 0, -E*A/L, 0, 0],\n", " [ 0, 12*E*I/L**3, 6*E*I/L**2, 0, -12*E*I/L**3, 6*E*I/L**2],\n", " [ 0, 6*E*I/L**2, 4*E*I/L, 0, -6*E*I/L**2, 2*E*I/L],\n", " [-E*A/L, 0, 0, E*A/L, 0, 0],\n", " [ 0, -12*E*I/L**3, -6*E*I/L**2, 0, 12*E*I/L**3, -6*E*I/L**2],\n", " [ 0, 6*E*I/L**2, 2*E*I/L, 0, -6*E*I/L**2, 4*E*I/L]])" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "def Mmat3D(L=1, A=1, rho=1):\n", " return rho*A*L/420 * np.matrix([\n", " [140, 0, 0, 0, 0, 0, 70, 0, 0, 0, 0, 0],\n", " [ 0, 156, 0, 0, 0, 22*L, 0, 54, 0, 0, 0, -13*L],\n", " [ 0, 0, 156, 0, -22*L, 0, 0, 0, 54, 0, 13*L, 0],\n", " [ 0, 0, 0, 140, 0, 0, 0, 0, 0, 70, 0, 0],\n", " [ 0, 0, -22*L, 0, 4*L**2, 0, 0, 0, -13*L, 0, -3*L**2, 0],\n", " [ 0, 22*L, 0, 0, 0, 4*L**2, 0, 13*L, 0, 0, 0, -3*L**2],\n", " [ 70, 0, 0, 0, 0, 0, 140, 0, 0, 0, 0, 0],\n", " [ 0, 54, 0, 0, 0, 13*L, 0, 156, 0, 0, 0, -22*L],\n", " [ 0, 0, 54, 0, -13*L, 0, 0, 0, 156, 0, 22*L, 0],\n", " [ 0, 0, 0, 70, 0, 0, 0, 0, 0, 140, 0, 0],\n", " [ 0, 0, 13*L, 0, -3*L**2, 0, 0, 0, 22*L, 0, 4*L**2, 0],\n", " [ 0, -13*L, 0, 0, 0, -3*L**2, 0, -22*L, 0, 0, 0, 4*L**2]])" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "def Kmat3D(L=1, A=1, E=1, Iy=1, Iz=1, GJ=1):\n", " a1 = E * A / L\n", " b1 = 12 * E * Iz / L**3\n", " b2 = 6 * E * Iz / L**2\n", " b3 = 4 * E * Iz / L\n", " b4 = 2 * E * Iz / L\n", " \n", " c1 = 12 * E * Iy / L**3\n", " c2 = 6 * E * Iy / L**2\n", " c3 = 4 * E * Iy / L\n", " c4 = 2 * E * Iy / L\n", " \n", " d1 = GJ / L\n", " \n", " return np.matrix([\n", " [ a1, 0, 0, 0, 0, 0, -a1, 0, 0, 0, 0, 0],\n", " [ 0, b1, 0, 0, 0, b2, 0, -b1, 0, 0, 0, b2],\n", " [ 0, 0, c1, 0, -c2, 0, 0, 0, -c1, 0, -c2, 0],\n", " [ 0, 0, 0, d1, 0, 0, 0, 0, 0, -d1, 0, 0],\n", " [ 0, 0, -c2, 0, c3, 0, 0, 0, c2, 0, c4, 0],\n", " [ 0, b2, 0, 0, 0, b3, 0, -b2, 0, 0, 0, b4],\n", " [ -a1, 0, 0, 0, 0, 0, a1, 0, 0, 0, 0, 0],\n", " [ 0, -b1, 0, 0, 0, -b2, 0, b1, 0, 0, 0, -b2],\n", " [ 0, 0, -c1, 0, c2, 0, 0, 0, c1, 0, c2, 0],\n", " [ 0, 0, 0, -d1, 0, 0, 0, 0, 0, d1, 0, 0],\n", " [ 0, 0, -c2, 0, c4, 0, 0, 0, c2, 0, c3, 0],\n", " [ 0, b2, 0, 0, 0, b4, 0, -b2, 0, 0, 0, b3]])\n" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "if el_type == pyaka._bernoulli_beam_2:\n", " kmat1 = Kmat2D(E=mat1.E, A=mat1.A, I=mat1.I, L=length)\n", " mmat1 = Mmat2D(A=mat1.A, L=length, rho=mat1.rho)\n", " nb_dofs = 3\n", "else:\n", " kmat1 = Kmat3D(E=mat1.E, A=mat1.A,\n", " Iy=mat1.Iy, Iz=mat1.Iz,\n", " GJ=mat1.GJ, L=length)\n", " mmat1 = Mmat3D(A=mat1.A, L=length, rho=mat1.rho)\n", " nb_dofs = 6" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "def assemble(Mg, Me, n1, n2, nd):\n", " s1 = n1 * nd\n", " s2 = n2 * nd\n", " \n", " Mg[s1:s1+nd, s1:s1+nd] += Me[ 0:nd , 0:nd ]\n", " Mg[s2:s2+nd, s2:s2+nd] += Me[nd:2*nd, nd:2*nd]\n", " Mg[s1:s1+nd, s2:s2+nd] += Me[ 0:nd , nd:2*nd]\n", " Mg[s2:s2+nd, s1:s1+nd] += Me[nd:2*nd, 0:nd ]" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "total_nb_dofs = nb_dofs * (nb_elem + 1)\n", "\n", "Ka = np.zeros((total_nb_dofs, total_nb_dofs))\n", "Ma = np.zeros((total_nb_dofs, total_nb_dofs))\n", "\n", "for e in range(nb_elem):\n", " n1, n2 = Conn[e, :]\n", " \n", " assemble(Ka, kmat1, n1, n2, nb_dofs)\n", " assemble(Ma, mmat1, n1, n2, nb_dofs)\n" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.581304729629375e-15\n", "4.999022389403736e-16\n" ] } ], "source": [ "print(np.linalg.norm(Ka - K.todense())/np.linalg.norm(Ka))\n", "print(np.linalg.norm(Ma - M.todense())/np.linalg.norm(Ma))" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "def analytical_solution(time, L, rho, E, A, I, F):\n", " def w(n):\n", " return n**2 * np.pi**2 / L**2 * np.sqrt(E * I / (rho * A));\n", "\n", " return 2 * F * L**3 / (np.pi**4 * E * I) * \\\n", " ((1. - np.cos(w(1) * time)) + (1. - np.cos(w(3) * time)) / 81. +\n", " (1. - np.cos(w(5) * time)) / 625.);\n" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "# Perform N time steps.\n", "# At each step records the displacement of all three nodes in x direction.\n", "N = 300\n", "\n", "disp = model.getDisplacement()\n", "disp[:, :] = 0.\n", "\n", "displs = np.zeros(N)\n", "\n", "for i in range(1, N):\n", " model.solveStep() \n", " displs[i] = disp[no_print, 1]\n" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEDCAYAAAA7jc+ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzlUlEQVR4nO3deVxU9f7H8deXXRRBBTfcMq3cMpXMTHPJyiXX1EpTU8vMJdtXq1u/9rrdMpfcl9TMzNBKy31fcQcUxQ0QWQQBUXa+vz+Y7iVDAR34zgyf5+Mxj4Y5Z855H869bw9nVVprhBBCOC4n0wGEEEKULCl6IYRwcFL0Qgjh4KTohRDCwUnRCyGEg5OiF0IIB2e06JVSc5RScUqpYCtNr45Sao1S6qhSKlQpVc8a0xVCCHtmeot+HtDVitNbAHyhtW4EtAbirDhtIYSwS0aLXmu9BUjM/5lS6lal1B9KqX1Kqa1KqTuKMi2lVGPARWu91jLtVK31FeunFkII+2J6i74gM4DxWutWwCvA1CJ+7zYgSSm1XCl1QCn1hVLKucRSCiGEnXAxHSA/pVQFoC3wk1Lqr4/dLcP6AR8U8LVzWuuHyVuW9kALIAL4EXgKmF2yqYUQwrbZVNGT9xdGktb6rqsHaK2XA8uv890o4IDW+hSAUioQaIMUvRCijLOpXTda6xTgtFJqAIDK07yIX98LVFJK+Vl+7gyElkBMIYSwK6ZPr/wB2AncrpSKUkqNBAYDI5VSh4AQoHdRpqW1ziFvn/56pdQRQAEzSya5EELYDyW3KRZCCMdmU7tuhBBCWJ+xg7G+vr66Xr16pmYvhBB2ad++fRe01n6Fj/k/xoq+Xr16BAUFmZq9EELYJaXU2eJ+R3bdCCGEg5OiF0IIBydFL4QQDk6KXgghHJwUvRBCODgpeiGEcHBS9EII4eBs7e6VQlxX8pUstp+8wKn4VADq+ZanfQM/vD1dDScTwnZJ0Qu7EBZziZlbT7HyYDSZObl/G+bm4kT3ptV5un19mvp7G0oohO2Sohc2LTLxCl+uCWPFwWjKuTrz2N216dPCnyY1KwIQej6FwAPn+OXAOQIPRtPjzhq8/OBt1PerYDi5ELbD2N0rAwICtNwCQVxLakY2k9afYN72MygFI++rxyin8/hs3wS7dkFMDFy8CElJ4OJCSuM7mdlxMLNULTJzYUiburz00G1U9JBdOsKxKKX2aa0DivUdKXpha7Ycj+fN5UeITk7j0btq8lLiAWpO+gKOHQMnJ7jzTqhbFypXBh8fyMiAQ4dgxw7ivX35ZsT7LHKtg5+XOx/2acpDTaqbXiQhrOZGil523QibobXm63Un+Gb9Cer7leentl4ETHwKDh+Gli1h4ULo1Qu8vAqeQEQEfh99xIeTnmdg/ea8+cS7jPp+H8Pvq8eb3Rrh5iInmYmySf6XL2xCelYOE5Yc5Jv1J+jf0p9VyZsI6N0REhNh2TIICoLBg69d8gB16sD06XD4MHdWdmX5h48yPCmEudvP8MTMXVxIzSi15RHClkjRC+MSUjN4ctZuVh6K5tX2tfhi9ut4/N/7MHQoHD0Kjz4KShV9go0awbZtuH/8Ee/Nfpsp22YSHHmRPlO2cyL2UsktiBA2SopeGBUel0rfqTs4fC6ZyQHlGTumF2rTJpgxA+bOhQo3ePaMszO89hps306P80dYuuBlMpJS6Dd1B1uOx1t1GYSwdVL0wphNYXH0m7qdyxnZLHEJ5ZHHu+QdbN2+HZ55pnhb8dfSujUcOEDz+1sSOGUU/hejGT53Dwt3FfvZDULYLSl6UepycjX/XhPGU3P3UtPThcAtk2j5/ivw2GNw8CAEFOuEgsJVrAiLFuH/zWcsW/AKHSIPMzEwmA9+DSUn18xZZ0KUJjnrRpSq+EsZTFhygB0nExjom80Hnz6BR0YazJ8PQ4ZYZyu+IErB8OFUuPtuZvbqxYcNH2YOjxAen8p/BjanSgX3kpmvEDZAtuhFqdl9KoEek7ay7+xFPk/dz+ev9sGjfj04cCDvwGtJlXx+TZvivGsX76Ue5uM/vmXX8Vi6T9rK7lMJJT9vIQyRohel4vudZxg0azflnTSBm75m4NT34I038vbHN2hQumGqVoX16xnU1Jdf5k7AM/ECT8zcxbRNJzF1AaEQJUmKXpSo3FzNJ6uO8s6KEDpVUaz8zzAaheyB33+HTz4BV0O3KPDwgO+/p8mYYfz67Qi6XzjGZ38c472VIbLfXjgc2UcvSozWmokrglm8O4IhVTL411uDcG7YAFasgFtvNR0vb1fRe+9R4dZb+XbkSGo+Mp4ZdCI5LYuvBt6Fs1Mp7EoSohRI0YsSobXm09XHWLw7gue8knjttSdRHTpAYGDe/WlsyZNPourW5a0+ffBOusAXDKCihysf9G6CKo3jBkKUMCl6USJmbzvN9C2neNIzmdcmPonq0weWLAF3Gz27pX172LmTsd27k6JcmU4fqlRw44Uut5lOJsRNK3QfvVLKQym1Ryl1SCkVopR6v4BxlFJqklIqXCl1WCnVsmTiCnuw8VgcH606SjePVD5470lU376wdKntlvxfbrsNdu7kjbRQ+h9Zx9frTvDroWjTqYS4aUU5GJsBdNZaNwfuAroqpdpcNU43oKHlNQqYZs2Qwn6ciL3E+B8O0Mglk39/PAynfn3hxx/NHXQtLj8/1Nq1fJwbRkBUCK/+sI/gc8mmUwlxUwotep0n1fKjq+V19WkJvYEFlnF3AT5KqRrWjSps3cXLmTy9IAiPnExmThqNZ7eH83bX2EvJ/8XTE7fAX5iWeYjKKRcYNXUT8ZfkzpfCfhXp9EqllLNS6iAQB6zVWu++ahR/IDLfz1GWz66eziilVJBSKig+Xm4s5UhyczUvLj3I+cQrTP/+Lfyb326fJf8Xd3f8Fs1jRvwWEtNzGDt5HdlXPatWCHtRpKLXWudore8CagGtlVJNrxqloFMT/nEystZ6htY6QGsd4OfnV+ywwnZN33KKTWHxvLNzIa08s2HlSihXznSsm+PqStN5k/k0dAV7kuGbX/aZTiTEDSnWBVNa6yRgE9D1qkFRQO18P9cC5ChWGbHvbCJfrgmjR/JJnty5PG9L3tZOobxR5cvT59t3GBC6kcl7Y9kRLn+JCvtTlLNu/JRSPpb35YAuwLGrRlsJDLWcfdMGSNZan7d2WGF70rNyeOWnw9RwyuKTOW+iPvsMmjc3Hcu6GjTg/Ydu5ZbEKF6cu4OU9CzTiYQolqJs0dcANiqlDgN7ydtH/5tSarRSarRlnFXAKSAcmAmMKZG0wuZMWn+C0xcu8+kvn1Pxvnvg+edNRyoRnmNH85/47cRlO/Hvn2UXjrAvhV4wpbU+DLQo4PPv8r3XwFjrRhO2LjQ6helbTjHg8inandgLy4/kPTjEESlF86//j6HPfckCutLv/iSa1/YxnUqIInHQ/1eKkqa15l+/huDjonl71lt5j+27/XbTsUpWvXq8fH8d/FITeWveNrn5mbAbUvTihmwMi2PP6URe2B+IT3VfeOst05FKRcWXX+Dt438Sclnx657TpuMIUSRS9KLYcnI1n60Oo55rNo+vnguffw6enqZjlQ43N3pOHEXj2JN8FbifzGw5t17YPil6UWzL90cRFnuJV9fMxLXNPdC/v+lIpcqpQwdedTtPBOVYukoOzArbJ0UviiU7J5dvN4TTTKXSfc/v8NVXpfMIQBvT8V/PExB9jElbzpCelWM6jhDXJUUviuXXw9FEJF5h/G/TUAMHwj33mI5khKpRg5frOxHn4snSHzebjiPEdUnRiyLLzdVM3XiSO3QqXUK3wfv/uGN1mdLmlWcIiAvnu6AYMmWrXtgwKXpRZGtCYzgRl8qY1TNwGjoE7rjDdCSjVIUKjG3hS7SHN4GzVpiOI8Q1SdGLItFa8+2GcG7JSaVH2HZ47z3TkWxCx+ceo2nyOaYGJ5Odlm46jhAFkqIXRbLpeDwh0Sk8t3YOziNHQL16piPZBOXiwri2tTjjVZXfv1tmOo4QBZKiF4XSWjN5Qzj+OZfpe2I7TJxoOpJNeWhQV25LjWPK8TRyMzJNxxHiH6ToRaF2nUpk39mLjN6wANfnRkPNmqYj2RQnZyfGtvTluHdN1sz42XQcIf5Bil4UasrGcPyyrzDg5A54/XXTcWxSj6HdqXslgSnByejsbNNxhPgbKXpxXQciLrIt/ALPbP0Bj/FjQZ4MViAXF2fGNKnIkUq12Tx7uek4QvyNFL24rikbw/HJTmfwqR3w8sum49i0vsMfoWZaEpOD4tA5cl69sB1S9OKaQqNTWHc0jhE7l1F+3HOO83jAEuLm5sKzDTwIqnILuxfIefXCdkjRi2uasikcr9xMhh3bAOPHm45jFx57pie+6SlM3hEFWu5XL2yDFL0o0Mn4VFYdPs+QPSvwHjkMKlc2HckueHi48UxNzbYqt3Lg5zWm4wgBSNGLa5i68STuOpuRh1bBSy+ZjmNXBo/ujU9GKlPWhpmOIgQgRS8KEJl4hcADUQzav4oqgwZA9eqmI9mVChXLM6JSGusq3Uromu2m4whReNErpWorpTYqpY4qpUKUUhMKGKejUipZKXXQ8nq3ZOKK0vDd5pM45+Yyau8v8OqrpuPYpWGje1EhM40pgfJgEmGeSxHGyQZe1lrvV0p5AfuUUmu11qFXjbdVa/2I9SOK0hSbks5PQZH0D15H9d7d5J42N8i7amWGlrvINNdbCN9xgAZtW5iOJMqwQrfotdbntdb7Le8vAUcB/5IOJsyYueUUOTm5PLdjKbzxhuk4dm3k6J64Z2cxbYnsvhFmFWsfvVKqHtAC2F3A4HuVUoeUUquVUk2sEU6UrsTLmSzafZbeYduo3bktNGpkOpJdq1K7OoOcYgn0qEPk4eOm44gyrMhFr5SqAPwMvKC1Trlq8H6grta6OfAtEHiNaYxSSgUppYLi4+NvMLIoKXO2nSY9K4cxWxfBW2+ZjuMQRo3sinNuLtPmrzcdRZRhRSp6pZQreSW/SGv9jxt5aK1TtNaplverAFellG8B483QWgdorQP85J4pNiU5LYv5O07T7XQQDQKaQKtWpiM5hOp33EL/rCiWqRrEnIw0HUeUUUU560YBs4GjWuuvrjFOdct4KKVaW6abYM2gomR9v/MMlzJyGLtpAbz9tuk4DuW5YZ3JcXJixoxVpqOIMqooW/T3AUOAzvlOn+yulBqtlBptGac/EKyUOgRMAh7XWq7/thdXMrOZve00nc8doclt/tC+velIDqV2y8b0STvL4mw/LpyNNh1HlEGFnl6ptd4GqELGmQxMtlYoUboW747g4pUsxm6YD7MK/KNN3KQxg+9n+S8RzJm2ktc+HV34F4SwIrkytoxLy8xh+uaT3Bt3nFbVykHXrqYjOaRb2zSn+5UIFmRWITk6znQcUcZI0Zdx83eeIT41k5fWzIJ//QvUdf94Ezdh3MB7SXXzZN6UX0xHEWWMFH0ZlpKexbRN4XQ6d4S7/StAz56mIzm0Ru1b0iU1gjmXfUiNvWA6jihDpOjLsFlbT5Ocls3La2bCRx/J1nwpGN+3FckeFZg/WR43KEqPFH0ZlZCaweytp+hxei9NG9eBBx4wHalMaP5Aa7pcOsv0VB+S4xJNxxFlhBR9GTVt00nSMrN5cd1s2ZovZS/2a0mKe3lmTwk0HUWUEVL0ZdD55DQW7DxDv7CtNGhzJ7RtazpSmdLkgTZ0Tw5nzqWKJMZdNB1HlAFS9GXQtxvC0Tk5TNgwDz780HScMunFvq247OrO9KkrTUcRZYAUfRlzNuEyS/dG8sSRtdR+6H5oIfdJN6HhQ+3oc/E481O9iJOtelHCpOjLmK/XncAlN5txmxfCu/IgMJMm9L+bLCcXpk773XQU4eCk6MuQsJhLBB44x7CDq6j6YAdo1sx0pDKt3kP30z8hlMWpFYiOka16UXKk6MuQr9aGUYEcRm9eDBMnmo4jgPED7kEDk2f8YTqKcGBS9GVESHQyf4bEMnL/Sip1aif3m7cRtbp24on4IyxNLU9EtGzVi5IhRV9GTN4QjpfKYfiWJfDOO6bjiHzGDrgX55xsvpm91nQU4aCk6MuAsJhLrA6OYfjBVXjfdw/ce6/pSCKfat0fYEjMfn65VI7w80mm4wgHJEVfBkzeGE55lcuITYtka94WKcXofvfgkZ3B13M3mE4jHJAUvYM7m3CZ3w5HMzR4DT6tmkOHDqYjiQL49u3B8Ihd/JbiztFzSabjCAcjRe/g5m4/g4vWDN9o2ZqXe9rYJicnRvUJoHzGFWYs3mI6jXAwUvQOLDkti6VBkfQ8s4eqjW6Fhx4yHUlch/cTAxgQuZdf4yEmKc10HOFApOgd2JI9EVzJzGHkxoWyNW8PnJ0Z0aUxOUqxYPFG02mEA5Gid1BZObnM23GGe+NP0KRGRXl6lJ2oM2IQD0ceZPHJK1zJzDYdRziIQoteKVVbKbVRKXVUKRWilJpQwDhKKTVJKRWulDqslGpZMnFFUa0OjuF8cjpPb14Mb7whW/P2ws2NkS2qkeRajp+XbjadRjiIomzRZwMva60bAW2AsUqpxleN0w1oaHmNAqZZNaUoFq01s7aeon5aAp2y46F/f9ORRDEEjBlM8/hTzNkXQ26uNh1HOIBCi15rfV5rvd/y/hJwFPC/arTewAKdZxfgo5SqYfW0okiCzl7kcFQyw7cswemFCeDiYjqSKAZVvjwj6zhz2t2Hjat2mo4jHECx9tErpeoBLYDdVw3yByLz/RzFP/8xQCk1SikVpJQKio+PL2ZUUVSztp7CJyedRyP3wciRpuOIG9Bt/BPUSE1g1ppQ01GEAyhy0SulKgA/Ay9orVOuHlzAV/7xN6fWeobWOkBrHeDn51e8pKJIziZcZk1ILIP3rMRz5FPg5WU6krgBrpUr8VTFVHZ61iBk+0HTcYSdK1LRK6VcySv5RVrr5QWMEgXUzvdzLSD65uOJ4sq7QCqXoYdWw/PPm44jbsLjY/rhmZnO3B+3mY4i7FxRzrpRwGzgqNb6q2uMthIYajn7pg2QrLU+b8WcoghS0rNYujeCnmHbqNarK/j/Y++ZsCPetarTzzmelW7+JBwNNx1H2LGibNHfBwwBOiulDlpe3ZVSo5VSoy3jrAJOAeHATGBMycQV1/PzviiuZOUyfNfPMH686TjCCoYNfZBMFzeWzPrNdBRhxwo9HUNrvY2C98HnH0cDY60VShRfbq7m+51nuSspkmbVK0BAgOlIwgoaNruVtmkbWZRdiWcvX8alfHnTkYQdkitjHcT2kxc4deEyQ7cthWeekQukHMiwtvWI9vJl7ZyVpqMIOyVF7yAW7DxLldwMup/dB4MHm44jrKhL/874X0lkXog8alDcGCl6B3AuKY31R2N57MBqPPr1AR8f05GEFTk7OzGkai67fepydNNe03GEHZKidwCLdp0FrRm8e0XebhvhcB4b3hX37AwWrthjOoqwQ1L0di4jO4cleyN54MJx/Gv5wn33mY4kSkAl/2p0z4phJVVJS0wyHUfYGSl6O7c2NJbEy5k8uWERPP20HIR1YAO7NOOSuyer5/xqOoqwM1L0dm7Zvihq5KbTLjoUhg41HUeUoDY97qPu5QR+PJ5sOoqwM1L0diw2JZ0tx+N59OCfOPfpDb6+piOJEqScnBhYQ7G78i2c2Sz76kXRSdHbseX7z5Gr4dGg3+UgbBnR/8mHcMrNYekvO0xHEXZEit5Oaa1Zti+Su1OiuKWSB3TqZDqSKAXV6lSjU2Ysy3J8yU65ZDqOsBNS9HbqQGQSJ+Mv03/7z3kHYZ1kVZYVA9s1IK58JTbNW2E6irAT0g52atm+KDx0Dt1P7ISnnjIdR5Sizn3uxzf9EksPx5qOIuyEFL0dSs/K4ddD0XQ/uQuvh7tADXlqY1ni6uJM30pZbKjUgIS9B03HEXZAit4O/RkSw6X0bPrvlYOwZdWjAzuQ7ezCyh/WmY4i7IAUvR1ati8K/4wU2ugkePhh03GEAXc0qkOzjASWJZeDtDTTcYSNk6K3M9FJaWw7cYFHg37DaeQIcHY2HUkY0r95dUL86hG68BfTUYSNk6K3Mz/vi0IDA0I2wogRpuMIg3r1vx/X3Gx+3nrcdBRh46To7UhuruanoEjanj9K7XtbQO3ahX9JOKxKFTzo4n6ZQO+GZIUeNR1H2DApejuy63QCERfTGChXwgqL/t1bkVDeh01zAk1HETZMit6O/BQUhVdOBl1Tz8Ajj5iOI2zA/a0b4pt9hWVRmZCRYTqOsFFS9HYiJT2LVYej6XV4PR4jngKXQp/rLsoAV2cn+t5SnvW17yLhp0DTcYSNKrTolVJzlFJxSqngawzvqJRKVkodtLzetX5MsfJgNBk5mseC1+fd8kAIi0f73ke2swsrVgWZjiJsVFG26OcBXQsZZ6vW+i7L64ObjyWu9tOeCO5IjKRZmybg7286jrAhd9T05k6nyywpdws6PNx0HGGDCi16rfUWILEUsohrOBaTwqHoFAYcWI0aPdp0HGGDnux0B8f96rJn1k+mowgbZK199PcqpQ4ppVYrpZpcaySl1CilVJBSKig+Pt5Ks3Z83+88i1tuNn2vnIUHHjAdR9ignh2aUDEngwWn0yEry3QcYWOsUfT7gbpa6+bAt0DgtUbUWs/QWgdorQP8/PysMGvHl5yWxfKgSHoHb6TyiCFyO2JRoHJuzgyo7cafdVoS97M8U1b83U23htY6RWudanm/CnBVSskz7azk531RpOVohoWsheHDTccRNuzJge3JdnZhyeoDpqMIG3PTRa+Uqq6UUpb3rS3TTLjZ6Yq8K2G/336KltFhNO3aDipXNh1J2LBbqnvTXiWxuEIDsk+fMR1H2JCinF75A7ATuF0pFaWUGqmUGq2U+uuoYH8gWCl1CJgEPK611iUXuezYGn6B0xfTGRa0EsaPNx1H2IEhD99JjJcv62bLjc7E/xR61Y3W+olChk8GJlstkfivBdtP4ZuWTDd/N2ja1HQcYQc6t29CzZUhzD+fQ9fMTHBzMx1J2AA5smejIhKusCEsnkH7f8ftedmaF0Xj4uzEkDsqsrNGI0IX/Gw6jrARUvQ2at6OMzjn5jIo8Sh07246jrAjgwY/QLnsDOZsPWU6irARUvQ26FJ6Fkt3n+GRo1uo/sxQebiIKBbv8u7098lgpV8T4jbtMB1H2AApehv0495IUrM1I0PWyCmV4oYMf7IzWc7OLPxxi+kowgZI0duY7Jxc5m49SeuoEJr16ADe3qYjCTtUv141HshNYKFrHdIjokzHEYZJ0duYNaGxnEvJZMSeQBg3znQcYcdG9GpJoqc3gdOXm44iDJOitzGzN4dTJyWOBxtVhdtvNx1H2LF72zWjUfoF5sS7o1NTTccRBknR25CDkUnsi0ph+J5fcH7lZdNxhJ1TSjGyTR2OV67F1qmLTccRBknR25DZW0/ilZXGAPckaNfOdBzhAHr2vx/fzFRmBSdBZqbpOMIQKXobEZ2UxqrD53n8wGoqvPg85N0+SIib4u7izNDbvNhSswkn5v5oOo4wRIreRszZdhq0ZljcQejb13Qc4UAGD+6Me04Wszceh9xc03GEAVL0NuDi5UwW7TxN75BN1Bo9XC6QElZVxcuDR6tqltdqRcxPK0zHEQZI0duAuTvOkJYDzx1fLxdIiRIxetgD5Dg5MWPFPpCby5Y5UvSGpWZkM2/zCR4O20HDscOhfHnTkYQDqlO1Ir19MllcowUJP680HUeUMil6w+ZvP01KNow5vQWefdZ0HOHAxgzvQoaLG7N+3Cb76ssYKXqDkq5k8t36MLqc2E3z8U+Bh4fpSMKBNajhQ49KOSyo1ZrExT+ZjiNKkRS9Qd9tCic1W/PK2c0wbJjpOKIMmDCsE1fcPJi+fA9kZ5uOI0qJFL0hsSnpzN16ij4hm7jj1THgUujDvoS4aQ1reNPbDxbUa0v8vEWm44hSIkVvyKS1YeTk5PJiwn4YMMB0HFGGPD+0Ixmubkz7/bBcLVtGSNEbcObCZX7cG8mgA6upM/EVcJLVIEpP/ape9KvhzMJb2xE7Y57pOKIUFNowSqk5Sqk4pVTwNYYrpdQkpVS4UuqwUqql9WM6lq/+PIprdibjMk7IYwKFEROGdCDXyZkp645DWprpOKKEFWVTch7Q9TrDuwENLa9RwLSbj+W4QqNTWHkklhF7Aqn6zutyTxthRO0q5RlQ150lDdpxbups03FECSu06LXWW4DE64zSG1ig8+wCfJRSNawV0NF8veYYFTMuM0qdgy5dTMcRZdi4wfeDk2Ly1rMg96t3aNbYOewPROb7OcrymbjKsZgU1hyLZ/jeFXi/97ZszQuj/H3K8fit5fmpYTsivpluOo4oQdYo+oLaqsCbaSilRimlgpRSQfHx8VaYtX2ZvDaMCplpDHe/AJ06mY4jBGOfaIezgklBcZCUZDqOKCHWKPoooHa+n2sB0QWNqLWeobUO0FoH+Pn5WWHW9uNkfCq/h8QyZN9v+Lz7lmzNC5tQraIHT97hzfKG93Hq33J4zVFZo+hXAkMtZ9+0AZK11uetMF2HMnVdGO7ZmYx0j4cOHUzHEeK/Rg+4F3dymRSSAhcumI7j8ILPJZOVU7r3GirK6ZU/ADuB25VSUUqpkUqp0Uqp0ZZRVgGngHBgJjCmxNLaqYiEKwQeOs/gA6vxfecN03GE+Bs/L3eGNvNlRcO2nPhiiuk4Di35ShaPT9/J+7+GlOp8C73uXmv9RCHDNTDWaokc0LT1x3DOzmaUezy0b286jhD/8Gy/1iw8soqvT2YyJSYGqlc3HckhLdh2ktTMHAbt+x36NCu1+colmSUsOimNZfujeezwn1Sb+JrpOEIUqHJ5N4a3rM7vDdsS+sm3puM4pCuZ2czZfILO4XtofN9dpTpvKfoSNmN9GDo3h2fdL0DbtqbjCHFNz/RqhVduFv855wwREabjOJwfdp/lYo4TYy8cgIcfLtV5S9GXoLhL6fywN5J+wRuoNfEV03GEuC5vT1eevseftQ3uIfSTSabjOJSM7Bxmrj3KPRFHaDXq8VI/606KvgTNWneMrFwY4x4H99xjOo4QhXqq211U0NlMOe8K4eGm4ziMX/afIyZTMfbUZhg4sNTnL0VfQhIvZ7JwdwS9jm6m3lsvmY4jRJF4e7oy5G5/Vt3elpMf/cd0HIeQnZPLtD+CufP8cdo/+YiRZ09I0ZeQOeuPkaYVY93j4O67TccRoshGdm2GO5qpce5w7JjpOHZvVXAMZ69oxoT8gRox3EgGKfoSkJyWxfwdZ+gWtoOGb71gOo4QxeJbwZ0nWtUksElHIj/+ynQcu6a1ZurqYBpciOCh3u3B09NIDin6EjB/YxiXcGacWwy0lNvzC/sz6uEmODkppse6wtGjpuPYrQ3H4jiWlMWYAytxGmvuWlIpeitLzchmztaTdDmxm8avyXVkwj7V8C5H/zurs/TOB4n9+EvTceyS1popf4RQKzmWnh0aQ+XKxrJI0VvZoi0nSNIujNNnISDAdBwhbthzDzchx9mFmdEKQkNNx7E7O08lsD82jWeDAnF96UWjWaTorSgtM4eZG47T/vR+7nrlWdNxhLgpdap40quxL4vu6kbiR5+bjmN3vv0jlKqpiQxo6ge1ahnNIkVvRUt2nOSCdmF8+gm5ClY4hDEPNyHd1Z05kbmyVV8MQWcS2Rl5iVF7luPx+qum40jRW0tGdg7T1+Rd+db6xRGm4whhFQ2redG1YWXmt+pJ8idfmI5jN779I5TKaSkMus0L7rjDdBwpemtZtusMMbkujE8JlvvNC4cyrlsTLrl7Mu90Bpw4YTqOzTsclcTmM8k8vecXPCe+ZToOIEVvFVmWK99anDvGfeOHyNOjhENpUtObB+t7MzugNymfyL76wnz7Ryje6akMqe9hE1vzIEVvFYFBEURluzA+4QDqoYdMxxHC6p7v3owUjwrMP54Kp0+bjmOzjp5PYW34RYYHrcBr4pum4/yXFP1NysnVTP3tEE1iwun0XOnflU6I0tCsljcP1KvIrIDeXPpUzqu/lil/hFAh8wrD67pCo0am4/yXFP1N+u1AFKezXBgfsxf1SA/TcYQoMRMeaUayhxcLQhIhMtJ0HJsTHpfK72EJDN33G95vv246zt9I0d+E3FzNlJX7uS3+LA8921+25oVDu7OWD53qVGBmq16kfv5v03FsztQ/QnDPymSkP9C4sek4fyNFfxPWBEdzPMOFsdG7cOrT23QcIUrchJ7NSSpXke8PxMD586bj2IyIhCusCIlj8KHVVJloW1vzIEV/w7TWfLs8iFsSz/HI031ka16UCXfV9qFDrfLMbNmLy1/InS3/Mu33QzjnZDOqjrPNbc1DEYteKdVVKRWmlApXSr1RwPCOSqlkpdRBy+td60e1LZuOxhKS7sJzUTtx7tvHdBwhSs3zPZuT6OnNwj2REBNjOo5x55LSWBZygceC11Pt3X/Uo00otOiVUs7AFKAb0Bh4QilV0D9ZW7XWd1leH1g5p03RWvP1sj34J8fS96ke4CR/GImyo1XdSrSv6cmMlj258v6HpuMYN2PlfnRuLqNvKwf16pmOU6CiNFRrIFxrfUprnQksAcr0Dun1R6I5dMWZ589uxbV/P9NxhCh1E3o1J8HTh0UHzpfp+9XHpqTzQ0gCjx7bjP/br5iOc01FKXp/IP+5VFGWz652r1LqkFJqtVKqSUETUkqNUkoFKaWC4uPjbyCuebm5mn8v3U29xGj6PddftuZFmRRQrzL31anI9NaPkvbG26bjGPPNkp3o3FzGNvOBatVMx7mmorRUQUcZ9VU/7wfqaq2bA98CgQVNSGs9Q2sdoLUO8PPzK1ZQW7F6dzhHs915IfkQro90Nx1HCGMmdGvCBU9vFkfnwObNpuOUupNxl/jxZCqDwjZT5/UJpuNcV1GKPgqone/nWkB0/hG01ila61TL+1WAq1LK12opbURGdg5fBh6k4YUIer7xtJxpI8q01rdU5t56PnzX9jHSXnsTcnNNRypVX85ej0dmOuO7NwNvb9NxrqsoRb8XaKiUukUp5QY8DqzMP4JSqrpSea2nlGptmW6CtcOaNuvH7ZxWnkx0j8K5xV2m4whh3EtdGxFfzpsp7rfC/Pmm45SaPUfPsTrZlWeiduH7zDDTcQpVaNFrrbOBccCfwFFgqdY6RCk1Wik12jJafyBYKXUImAQ8rrW+eveOXYuKv8S3BxPodiaIDh/a7kEXIUrT3fUq069FTaa3GUD4R19BgsNt3/1DelYOb87dhn9yLM+M7wfOzqYjFUqZ6uOAgAAdFBRkZN7FpbXmmXeXsP2yG+tagf/gR01HEsJmJKRm8MDnG7jt5GF+9AxHzZxpOlKJ+mrmGiadzGJe5n46fvVOqc9fKbVPa12sB1LLKSNFsOSnrazLqsjLSYfwHySnUwqRX5UK7rzZswl7ajdl7oE42LTJdKQSExp+nmnH0+gTEUTHD182HafIpOgLER5+jvf3JNAu5igjPn9eDsAKUYCBAbXp0rAKn3YaQfALb0NysulIVpealsm4qRuolJbMO093Bk9P05GKTIr+OtLSsxg/eT3lstL597C2OPk63IlEQliFUorPH29J5fJujG89jNTxL5iOZFVaayb+3w+ccfFiUp00qnRubzpSsUjRX4PWmtc/+pFj7pX4quYlqnVqazqSEDatcnk3vh7amrOVa/JSqj+5S340HclqFn63gsBcXyZcOUabV0eZjlNsUvTXMGPWH6zMqsQrSYfo9NozpuMIYRfa1K/CxO6NWHPbvXw9ew1ERJiOdNN2/LGLf512onPCccZ9Pt4ud99K0Rdg89ZgPjuRTY/ow4z59wt2uWKFMGV4+/oMuK0ik1r1Y/Uzb0JqqulIN+zMgWM892cE9VNi+eat/jh7ljMd6YZI0V8lLPw84wKPcVtiJF+80Q9VoYLpSELYFaUUHw5tS8uK8FKTfgSPnGCXV82mREbz9MztKJ3L7BH34HVL7cK/ZKOk6POJjrnIsCmb8Uy/wqwe9fBsfLvpSELYJXcXZ74b/wCVyzkzrGpnTj/zPGRnm45VZGkXk3n6o0DOlPdlauca1Gl9p+lIN0WK3iI2JpEhn/7GZZyYF+BBrb5ywzIhbkZVLw++n/AAunx5hri2IHrI05CZaTpWoTIvJjHm9Xns9a7NV42caNurg+lIN02KHog8epoBH/9OjJMHs+qn02j4QNORhHAI9at6MW9cR5Ir+dG/UgfC+wyCixdNx7qm+KhYnnxzMRsrN+Cjetn0GtHTdCSrKPNFv3vtHvp9t5NkJzcW3l2Oe8YPNR1JCIdyZy0flozvQGZlX/o37M/uRwbDyZOmY/3D/iNneOTLDRwuX52vG+Yy6Lm+piNZTZkt+txczZRpv/PE2lgqZKWztE99WjzxiOlYQjikJjW9Wf7iA1T29Wbwfc8yf/jb6K1bTcf6r8Vrj/DY94dwy0jj5zbl6DPSMbbk/1Imiz4hNYOn3l/GF2ehR1wwv772ILd3vNt0LCEcWp0qngS+2oWO9bx5r+0QXv5iBekLFxvNlJ6Vwxvfreet9RHcGxnMr71q02RAN6OZSkKZK/rgiER6fPAbuy678PH5LUyaPJ4K9ez3tCkh7ElFD1dmjO7AC/fVYnnjTvRfF0/Ue5+Agbvonk9O47GPfmXJmXTGhfzB3Fe749PtwVLPURrKVNFv2H+Ggd9uwSk1heX6IIPmfCznyQtRypycFC/0bM6sQXdxtmpdel2sy46RL0F6eqllCI26SJ+PV3EyOYvp4St5Ze57OLdsUWrzL21lpugXrDrA0z8eoX78WX5pnEXTf79vFw8MEMJRdbnTn8BXHqByBXeG+HZi1tA30XFxJT7fLYfOMvCbTajUVJZxiIeXTAEHv2Ghwxd9Tq7m/2Zv5N0t0XQ+e4Clj91BtXH2d1MiIRzRrVW9CHyvNw9W1nxYvwuvjPmGzHXrS2x+S/88yIhFh6iVcI5fGmVyx38+BBeXEpufrXDook/NyOa5TwOZfeIKT4VtZPpbffHs9rDpWEKIfCq4uzDt9V680LQiPzdoy4hZO0l57U2r7srJydV8sWALr208x73nQvlpwG3UGFd2blbosEUfHnGB3hN/Zl2SM++eXs+/vnsV5zubmY4lhCiAUooXnmzPl73uYFfd5vRKqk9wxx6we/dNTzs5LYunv/qDKaGXeDx8G3Ne74FXt4eskNp+ONzfLDm5mgUL1/PF4RQ8M7JY6HGCtgs/A1dX09GEEIXo3/ZW6tasxPi5O+nnPZ5XJnzO8I634fqv98DDo9jTWx8aw9uL95CQkcuHIb8y+Lv3UHXqlEBy2+ZQW/RBR87S7/XFvB+aQUD8SX7rWp22X0yUkhfCjtxdrzKrXnuQ+++oxscdR9Azzp9dXR6FbduKPI2zCZcZM3cnIxfsw/t8JMuO/ciTi74skyUPoHQRzl9VSnUFvgGcgVla60+vGq4sw7sDV4CntNb7rzfNgIAAHRQUdKO5/yb8TCyfzdnI2kwv/FITmegWRa/3x6G8vKwyfSFE6dNa82dILB/8tI/oDOh4MohXvRJo8vpYuL3gO8smXs7k29XBLAw6h0tWJs/tWc7obs1we/1VhznLTim1T2sdUJzvFLrrRinlDEwBHgSigL1KqZVa69B8o3UDGlpe9wDTLP8tMdkZmexevYMFO06zFl88s5x5OXEnI59/FM9WQ0py1kKIUqCUomvT6nS8vSvzNx1nqs6ih3Kly9uLeMw5no5dWuF6f3u0hwdhR06x6MB5frrkSSZOPHZ4LS+Wj6fq5DeheXPTi2JcUfbRtwbCtdanAJRSS4DeQP6i7w0s0Hl/HuxSSvkopWporc9bO/DGhav4v52xnPeoSJqrB5XSPXg2M5SnH29HlY4DrD07IYRhHq7OPPtgIx5v14DZfwSzmBasww33YxlUD9pCkocXyeW8cMv2pM/ZvYyqrWjw8XBo2dJ0dJtRlKL3ByLz/RzFP7fWCxrHH/hb0SulRgGjAOrc4L4ynyoVaaRO05ErBNSvQee+D+HhV+WGpiWEsB/e5Vx5qW8LxvdqzqZjcezeE0Z8DFTQ2TSrBl06tMD31t7y6M8CFKXoC/qtXb1jvyjjoLWeAcyAvH30RZj3P7To1o4p3drdyFeFEA7A1dmJB5tU58Em1U1HsRtFOesmCsh/169aQPQNjCOEEMKAohT9XqChUuoWpZQb8Diw8qpxVgJDVZ42QHJJ7J8XQghRfIXuutFaZyulxgF/knd65RytdYhSarRl+HfAKvJOrQwn7/TK4SUXWQghRHEU6cpYrfUq8so8/2ff5XuvgbHWjSaEEMIaHOrKWCGEEP8kRS+EEA5Oil4IIRycFL0QQji4It3UrERmrFQ8cPYGv+4LXLBiHJNkWWyTLIttkmWBulprv+J8wVjR3wylVFBx795mq2RZbJMsi22SZbkxsutGCCEcnBS9EEI4OHst+hmmA1iRLIttkmWxTbIsN8Au99ELIYQoOnvdohdCCFFEUvRCCOHotNal8gK6AmHk3eHyjQKGK2CSZfhhoGVh3wUqA2uBE5b/Vso37E3L+GHAw/k+bwUcsQybhGX3lZ0uyybLZwctr6q2vCxAFWAjkApMvmo+drVeClkWe1svDwL7LL//fUBnO14v11sWe1svrfNlPQT0vdH1UqyFvNEXebc3PgnUB9wsoRtfNU53YLXlF9UG2F3Yd4HP//qFAW8An1neN7aM5w7cYvm+s2XYHuBey3xWA93seFk2AQF2tF7KA+2A0fyzHO1tvVxvWextvbQAalreNwXO2fF6ud6y2Nt68QRcLO9rAHH5fi7WeimtXTf/fcC41joT+OsB4/n99wHjWutdgI9SqkYh3+0NzLe8nw/0yff5Eq11htb6NHn/6rW2TK+i1nqnzvttLcj3HbtalmJmtoll0Vpf1lpvA9Lzz8Ae18u1lsVKSntZDmit/3oiXAjgoZRyt9P1UuCyFDOzrSzLFa11tuVzDyyPZ72R9VJaRX+th4cXZZzrfbeatjzJyvLfqkWYVlQhOexlWf4yVyl1UCn1jlLFfipyaS/L9XLY23opjL2ul0eBA1rrDOx/veRflr/Y1XpRSt2jlAohbzfNaEvxF3u9lFbR38wDxov04PESnFZRp11S87/edwZrrZsB7S2vIYVMqzjTLmwca/wui5PDGtMojWUBO10vSqkmwGfAs8XIUehkizCN0lgWsMP1orXerbVuAtwNvKmU8riRaZVW0d/MA8av991Yy58xf/05E1eEadUqJIe9LAta63OW/14CFlP8XTqlvSzXy2Fv6+Wa7HG9KKVqAb8AQ7XWJ/PNw+7WyzWWxS7XS77sR4HL5B13KP56uXqnfUm8yHtk4SnyDib+dSCiyVXj9ODvBzH2FPZd4Av+fhDjc8v7Jvz9AOYp/ncAc69l+n8dxOhuj8timZavZRxXYBl5f9rZ7LLkm+ZT/PMApl2tl2stiz2uF8DHMt6jBWSxq/VyrWWx0/VyC/87+FqXvDL/axmKtV5KpegtwboDx8k78vy25bPRf/2yLYGnWIYfId/R8YK+a/m8CrCevNOS1gOV8w172zJ+GPmOSAMBQLBl2GRu7HQx48tC3lkf+8g7hSsE+AbLP2Y2vixngETyTkuM4n9nHtjjevnHstjjegEmkre1eJCrTj20t/VyrWWx0/UyxJL1ILAf6HOjPSa3QBBCCAcnV8YKIYSDk6IXQggHJ0UvhBAOTopeCCEcnBS9EEI4OCl6IYRwcFL0Qgjh4P4fgCovdARvPdEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def sol(x):\n", " return analytical_solution(x, L, mat1.rho, mat1.E,\n", " mat1.A, mat1.I, F)\n", "\n", "times = np.arange(N) * deltaT\n", "plt.plot(times, displs, 'r')\n", "plt.plot(times, sol(times))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What I do not fully understand is why the middle node first go backwards until it goes forward.\n", "I could imagine that there is some vibration, because everything is in rest." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.max(displs - sol(times))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.1+" } }, "nbformat": 4, "nbformat_minor": 4 }