386 lines
38 KiB
Plaintext
386 lines
38 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "e465f900-0a94-46d4-b3c5-23cc0129557b",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy\n",
|
|
"\n",
|
|
"from ipywidgets import IntSlider, interact\n",
|
|
"from typing import Iterator, Tuple"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "b39df788-2971-4974-9ced-1518a47181e8",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"class ELM:\n",
|
|
" n_rows: int\n",
|
|
" n_cols: int\n",
|
|
" displacement: numpy.ndarray\n",
|
|
" velocity: numpy.ndarray\n",
|
|
" force: numpy.ndarray\n",
|
|
"\n",
|
|
" # K\n",
|
|
" elastic_spring_constant: float\n",
|
|
" # c\n",
|
|
" bond_bending_constant: float\n",
|
|
" \n",
|
|
" def make(n_rows: int, n_cols: int):\n",
|
|
" return ELM(\n",
|
|
" n_rows,\n",
|
|
" n_cols,\n",
|
|
" numpy.zeros((n_rows, n_cols, 2)),\n",
|
|
" numpy.zeros((n_rows, n_cols, 2)),\n",
|
|
" numpy.zeros((n_rows, n_cols, 2))\n",
|
|
" )\n",
|
|
" \n",
|
|
" def __init__(\n",
|
|
" self,\n",
|
|
" n_rows: int,\n",
|
|
" n_cols: int,\n",
|
|
" displacement: numpy.ndarray,\n",
|
|
" velocity: numpy.ndarray,\n",
|
|
" force: numpy.ndarray\n",
|
|
" ) -> None:\n",
|
|
" self.n_rows = n_rows\n",
|
|
" self.n_cols = n_cols\n",
|
|
" self.displacement = displacement\n",
|
|
" self.velocity = velocity\n",
|
|
" self.force = force\n",
|
|
"\n",
|
|
" self.elastic_spring_constant = 1.0\n",
|
|
" self.bond_bending_constant = 2.0\n",
|
|
" \n",
|
|
" # u_i(t)\n",
|
|
" def get_displacement(self, row: int, col: int) -> numpy.ndarray:\n",
|
|
" return self.displacement[row, col, :]\n",
|
|
"\n",
|
|
" def set_displacement(self, row: int, col: int, displacement: numpy.ndarray) -> None:\n",
|
|
" self.displacement[row, col, :] = displacement\n",
|
|
"\n",
|
|
" # x_i(t)\n",
|
|
" def get_lattice_position(self, row: int, col: int) -> numpy.ndarray:\n",
|
|
" return numpy.array([col, row], dtype=float)\n",
|
|
"\n",
|
|
" # r_i(t)\n",
|
|
" def get_actual_position(self, row: int, col: int) -> numpy.ndarray:\n",
|
|
" return self.get_lattice_position(row, col) + self.get_displacement(row, col)\n",
|
|
" \n",
|
|
" # v_i(t)\n",
|
|
" def get_velocity(self, row: int, col: int) -> numpy.ndarray:\n",
|
|
" return self.velocity[row, col, :]\n",
|
|
"\n",
|
|
" # F_i(t)\n",
|
|
" def get_force(self, row: int, col: int) -> numpy.ndarray:\n",
|
|
" return self.force[row, col, :]\n",
|
|
"\n",
|
|
" def set_force(self, row: int, col: int, force: numpy.ndarray) -> None:\n",
|
|
" self.force[row, col, :] = force"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "fb5633c4-c2c0-4a9c-bde8-8a99c5db8a9d",
|
|
"metadata": {},
|
|
"source": [
|
|
"(O'Brien, 2008)\n",
|
|
"\n",
|
|
"The force acting on node $i$ is given by\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\mathbf{F}_i = \\sum^N_j \\mathbf{F}_{ij} \\frac{\\mathbf{r}_{ij}}{|\\mathbf{r}_{ij}|}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Where $N = 8$ (for the 2D case) is the number of neighbors $j$.\n",
|
|
"The force $\\mathbf{F}_{ij}$ on node $i$ from node $j$ is given by\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\mathbf{F}_{ij} = -K_{ij} (\\mathbf{u}_{ij} \\cdot \\mathbf{x}_{ij}) + \\frac{c \\mathbf{u}_{ij}}{|\\mathbf{x}_{ij}|^2}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"where\n",
|
|
"\n",
|
|
"* $K_{ij} \\in \\mathbb{R}$ is the elastic spring constant\n",
|
|
"* $\\mathbf{u}_{ij} = \\mathbf{u}_i - \\mathbf{u}_j \\in \\mathbb{R}^2$ is the displacement\n",
|
|
"* $\\mathbf{x}_{ij} = \\mathbf{x}_i - \\mathbf{x}_j \\in \\mathbb{R}^2$ is the vector connecting nodes $i$ and $j$ on the undistorted lattice\n",
|
|
"* $c \\in \\mathbb{R}$ is the bond-bending constant\n",
|
|
"\n",
|
|
"The $\\frac{c \\mathbf{u}_{ij}}{|\\mathbf{x}_{ij}|^2}$ part does not make sense to me since the first summand in $\\mathbf{F}_{ij}$ is scalar so the implementation below leaves it out."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "b4ed52bd-02aa-411a-9af4-be9bb28227bc",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlyElEQVR4nO3de2yT9+Hv8Y9DiKEQOw2F2BkhTW9cikgOaRus3tZCm/Y3QRggtT+xiW3oVGUBNYH+tPJHm6Jfp3CKtK2s0Ms2lUkb0FKJVvS3duMECOoGlHJRoS0RdGhEShzWtbmQNQ6HPOcPFo9wCTaxv7af7/slRRr2O+bZU6Q8SvzJ43EcxxEAAIAhWak+AAAAYBcuPgAAgFFcfAAAAKO4+AAAAEZx8QEAAIzi4gMAABjFxQcAADCKiw8AAGBUdqoP4GJ9fX1qaWlRbm6uPB5Pqg8HAADEwHEcdXV1qbCwUFlZg39vI+0uPlpaWlRUVJTqwwAAANegublZ48ePH7RJu4uP3NxcSecP3ufzpfhoAABALDo7O1VUVBT9Oj6YtLv46P9Ri8/n4+IDAIAME8tbJnjDKQAAMIqLDwAAYFRcFx/PP/+8PB7PgI9JkyZFn+/p6VF1dbXGjBmj0aNHa/78+Wpra0v4QQMAgMwV93c+br/9drW2tkY/Pvzww+hztbW12rZtm7Zs2aLGxka1tLRo3rx5CT1gAACQ2eJ+w2l2drYCgcAlj3d0dOg3v/mNNm7cqAcffFCS9MYbb2jy5Mnau3evZsyYMfSjBQAAGS/u73wcP35chYWFuummm7Rw4UKdOnVKknTgwAGdPXtWs2bNiraTJk3ShAkTtGfPnsQdMQAAyGhxfeejoqJCGzZs0MSJE9Xa2qpVq1bp3nvv1dGjRxUOh5WTk6O8vLwBn1NQUKBwOHzF14xEIopEItE/d3Z2xvf/AAAAZJS4Lj4effTR6P+eNm2aKioqVFxcrLfeeksjR468pgOor6/XqlWrrulzAQBA5hnSLxnLy8vTbbfdphMnTuihhx5Sb2+v2tvbB3z3o62t7bLvEem3cuVKLV++PPrn/t+Qlirn+hx9dPIrne7q0bjcEbqrJF/Dsjw0NDQ0NDRWNCYM6eLjzJkz+uKLL/T9739f5eXlGj58uBoaGjR//nxJUlNTk06dOqVQKHTF1/B6vfJ6vUM5jIT54GirVm37TK0dPdHHgv4Rqps9RY9MDdLQ0NDQ0Li6McXjOI4Ta/z0009r9uzZKi4uVktLi+rq6nT48GF99tlnGjt2rJYsWaI//OEP2rBhg3w+n5YtWyZJ+stf/hLzAXV2dsrv96ujo8Por1f/4GirlvzuoC4+Gf3Xg698b7ok0dDQ0NDQuLIZ6gVIPF+/47r4ePzxx7V792794x//0NixY3XPPffopz/9qW6++WZJ53/J2IoVK7Rp0yZFIhFVVlZq/fr1g/7YZSgHnyjn+hzd8392DLgavJBHUoHPK8mjcCcNDQ0NDY27moB/hD78yYND+hFM0i4+TEjFxceeL/6h//zVXiN/FwAA6WjT/56h0M1jrvnz4/n6zb1dJJ3uuvzVIAAAtjD5tZCLD0njckek+hAAAEgpk18Lh7R2cYu7SvIV9I9QuKPnkjfjSAN/ZtbWSUNDQ0ND464m4D8/uzWF73xIGpblUd3sKZf9jyJJjqTn59yu5+fQ0NDQ0NC4r6mbPcXo7/vg4gMAABjF2kVMbWloaGho7G6Y2jK1BQDAOKa2hjG1BQDYjqmtYUxtAQC2Y2prGFNbGhoaGhqbG6a2KcDUloaGhobG5oapLQAAcDXWLmJqS0NDQ0Njd8PUlqktAADGMbU1jKktAMB2TG0NY2oLALAdU1vDmNrS0NDQ0NjcMLVNAaa2NDQ0NDQ2N0xtAQCAq7F2EVNbGhoaGhq7G6a2TG0BADCOqa1hTG0BALZjamsYU1sAgO2Y2hrG1JaGhoaGxuaGqW0KMLWloaGhobG5YWoLAABcjbWLmNrS0NDQ0NjdMLVlagsAgHFMbQ1jagsAsB1TW8OY2gIAbMfU1jCmtjQ0NDQ0NjdMbVOAqS0NDQ0Njc0NU1sAAOBqrF3E1JaGhoaGxu6GqS1TWwAAjGNqaxhTWwCA7ZjaGsbUFgBgO6a2hjG1paGhoaGxuWFqmwJMbWloaGhobG6Y2gIAAFdj7SKmtjQ0NDQ0djdMbZnaAgBgHFNbw5jaAgBsx9TWMKa2AADbmfxayMWH/j21HUzA51XAR0NDQ0ND474myNTWvGFZHs0pDQ7aVJUVqqqMhoaGhobGfc2c0qDRqS1vONXV1y7S+SvHwd4tTENDQ0NDk6lNkLULaxcAAExj7WIYaxcAgO1YuxjG2gUAYDtuLGcYN5ajoaGhobG54cZyKcCN5WhoaGhobG64sRwAAHA11i7ixnI0NDQ0NHY33FiOqS0AAMYxtTWMqS0AwHZMbQ1jagsAsB1TW8OY2tLQ0NDQ2Nxk1NR29erV8ng8qqmpiT7W09Oj6upqjRkzRqNHj9b8+fPV1tY21ONMKqa2NDQ0NDQ2Nxkztd2/f79ee+01TZs2bcDjtbW12rZtm7Zs2aLGxka1tLRo3rx5Qz5QAADgDte0djlz5oymT5+u9evX64UXXlBZWZl+8YtfqKOjQ2PHjtXGjRu1YMECSdKxY8c0efJk7dmzRzNmzLjqazO1paGhoaGhYWp7iUWLFik/P18///nP9e1vfzt68bFjxw7NnDlTX3/9tfLy8qJ9cXGxampqVFtbe8lrRSIRRSKRAQdfVFTE1BYAAINMTm3jfsPp5s2bdfDgQe3fv/+S58LhsHJycgZceEhSQUGBwuHwZV+vvr5eq1ativcwEoqpLQDAdmk7tW1ubtZTTz2l3//+9xoxIjGTnJUrV6qjoyP60dzcnJDXjQdTWwCA7Ux+LYzr4uPAgQM6ffq0pk+fruzsbGVnZ6uxsVFr165Vdna2CgoK1Nvbq/b29gGf19bWpkAgcNnX9Hq98vl8Az5M65/aXuknXR5JAZ9XAR8NDQ0NDY37mmA6T21nzpypI0eO6PDhw9GPO+64QwsXLoz+7+HDh6uhoSH6OU1NTTp16pRCoVDCDz5RmNrS0NDQ0NjcmJ7axvWej9zcXE2dOnXAY6NGjdKYMWOijy9evFjLly9Xfn6+fD6fli1bplAoFNPSBQAAuN+Qbyx34dpFOv9LxlasWKFNmzYpEomosrJS69evv+KPXS7G1JaGhoaGhoaprVHc1RYAAPO4q61hTG0BALZL26mtWzG1BQDYjrvaGsZdbWloaGhobG4y6q62bsHUloaGhobG5iZj7moLAABwLVi7iKktDQ0NDY3dDVNbprYAABjH1NYwprYAANsxtTWMqS0AwHZMbQ1jaktDQ0NDY3PD1DYFmNrS0NDQ0NjcMLUFAACuxtpFTG1paGhoaOxumNoytQUAwDimtoYxtQUA2I6prWFMbQEAtmNqaxhTWxoaGhoamxumtinA1JaGhoaGxuaGqS0AAHA11i5iaktDQ0NDY3fD1JapLQAAxjG1NYypLQDAdkxtDWNqCwCwncmvhVx86N9T28EEfF4FfDQ0NDQ0NO5rgkxtzRuW5dGc0uCgTVVZoarKaGhoaGho3NfMKQ0andryhlNdfe0inb9yHOzdwjQ0NDQ0NJnaBFm7sHYBAMA01i6GsXYBANiOtYthrF0AALbjxnKGcWM5GhoaGhqbG24slwLcWI6GhoaGxuaGG8sBAABXY+0ibixHQ0NDQ2N3w43lmNoCAGAcU1vDmNoCAGzH1NYwprYAANsxtTWMqS0NDQ0Njc0NU9sUYGpLQ0NDQ2Nzw9QWAAC4GmsXMbWloaGhobG7YWrL1BYAAOOY2hrG1BYAYDumtoYxtQUA2I6prWFMbWloaGhobG6Y2qYAU1saGhoaGpsbprYAAMDVWLuIqS0NDQ0Njd0NU1umtgAAGMfU1jCmtgAA2zG1NYypLQDAdkxtDWNqS0NDQ0Njc8PUNgWY2tLQ0NDQ2NwwtQUAAK7G2kVMbWloaGho7G6Y2jK1BQDAOKa2hjG1BQDYjqmtYUxtAQC2Y2prGFNbGhoaGhqbm7Se2r7yyiuaNm2afD6ffD6fQqGQ3n///ejzPT09qq6u1pgxYzR69GjNnz9fbW1tCT/oRGNqS0NDQ0Njc5PWU9vx48dr9erVOnDggD7++GM9+OCDqqqq0qeffipJqq2t1bZt27RlyxY1NjaqpaVF8+bNS8qBAwCAzDTktUt+fr7WrFmjBQsWaOzYsdq4caMWLFggSTp27JgmT56sPXv2aMaMGTG9HlNbGhoaGhoapraXde7cOW3ZskWLFi3SoUOHFA6HNXPmTH399dfKy8uLdsXFxaqpqVFtbe1lXycSiSgSiQw4+KKiIqa2AAAYlNZT2yNHjmj06NHyer168skntXXrVk2ZMkXhcFg5OTkDLjwkqaCgQOFw+IqvV19fL7/fH/0oKiqK95CGjKktAMB2aT21nThxog4fPqx9+/ZpyZIlWrRokT777LNrPoCVK1eqo6Mj+tHc3HzNr3WtmNoCAGyX1lPbnJwc3XLLLZKk8vJy7d+/Xy+99JIee+wx9fb2qr29fcB3P9ra2hQIBK74el6vV16vN/4jTyCmtjQ0NDQ0NjdpPbW9nL6+PkUiEZWXl2v48OFqaGiIPtfU1KRTp04pFAoN9a9JqmFZHs0pDV5xhiSl3yyKhoaGhoYmU6e2cX3nY+XKlXr00Uc1YcIEdXV1aePGjdq1a5f++Mc/yu/3a/HixVq+fLny8/Pl8/m0bNkyhUKhmJcuqfLB0Va9vvvkFZ9/4r4SPTI1qA+Otho8KgAA3CmutcvixYvV0NCg1tZW+f1+TZs2TT/5yU/00EMPSTr/S8ZWrFihTZs2KRKJqLKyUuvXrx/0xy4XMz21vdrMVpKC/hFq/K8HdP+ancxxaWhoaGhc12TM1DZZTF98xDqzffY7k/Xf//N50o8HAIBUSOuprdvEOi3621f/TPKRAACQOmk9tXWbWKdFxfnXJflIAABIHZNTW+svPr7ujmiwH3F5dP49H98P3aigf/D/MAGfVwEfDQ0NDQ1NZjXBTJvaZrIPjraqeuMh9V3lXS91s6coJztLc0qDg3ZVZYWqKqOhoaGhocmsZk5p0OjU1to3nMaycsnySC//53T9x7RgTH3gKu8opqGhoaGhSccmyNrFzMVHrCuX/nf/cvM5AICbsXYxINZ39fZ33HwOAOBmrF0MiPVdvf0dN58DALgZaxcDyouvV/6onCs+379y6X/3b//N56700zCP/v2OYhoaGhoamkxqWLsY8MHRVt2/Zqe+6u697PP9/4EuvNHOsCyP6manz02AaGhoaGhorLixnBt8cLRVS3538Ir/ESQp4B+hutlT9MjUwadJAAAgflatXWKZy+aPGq69K2cpJ3vgN4Wu9rnpeKMgGhoaGhoabiwXg2RefMQ7r72WzwUAIBMxtU2SeOe11/K5AABkIqa2SRLvvPZaPhcAgEzE1DZJ4p3XXoipLQ0NDQ2NWxumtklyLfPaCzG1paGhoaFxa8PUNgmY1wIAkD5cv3YZyrw2ntdJx+kUDQ0NDQ0NU9sYJPriYyjz2mt5HQAAMhFT2wQayrz2Wl4HAIBMxNQ2gYYyr72W1wEAIBMxtU2gWCaysUyMmNrS0NDQ0Li1YWqbYMOyPJpTGhx06RLLxIipLQ0NDQ2NWxvTU1vXX3x8cLRVr+8+ecXnn7ivhHktAAAGuXrtEsvMNhjjvIipLQ0NDQ2NWxumtgm8+EjUzDae1wIAIBMxtU2QRM1s43ktAAAyEVPbBEnUzDae1wIAIBMxtU2QRM1sY32tdJtO0dDQ0NDQMLU1rH8eeyWOpDmlwZjeYMPUloaGhobGrQ1T2wR7ZGpQT9xXcsXnX999Uh8cbTV4RAAA2M3VaxcptolsLBMjprY0NDQ0NG5tmNpyV1sAAIxjaptA3NUWAICrY2qbQNzVFgCAq2Nqm0BXm8hK0phROSovvn5Ir5OO0ykaGhoaGhqmtilw4dz2Sif+H929un/NzkFXL0xtaWhoaGjc2jC1TYJHpgb1yvemK+C/8reUwh09WvK7g8xuAQBIMtevXS7U+//6NKP+/+qr7rOXfX6wuRFTWxoaGhoatzZMbZN48TGU2S1TWwCAmzG1TZKhzG6Z2gIA3IypbZIMZXbL1BYA4GZMbZNkKLNbprY0NDQ0NG5tmNom0VBmt0xtaWhoaGjc2jC1TTJmtwAApJZVa5cLxTu7ZWpLQ0NDQ+PWhqmtoYuPeGe3TG0BAG7G1NaAeGe3TG0BAG7G1NaAeGe3TG0BAG7G1NaAWGa3WR7p6+7IgH4w/XMmGhoaGhqaTGqY2hpy4ez2SvocqXrjIX1wtFXDsjyaUxoctK8qK1RVGQ0NDQ0NTWY1c0qDRqe21r7htN8fPmnV0k0H1XeFs9D/LuDG/3pA96/ZecW1i3T+6nKwdxTT0NDQ0NCkYxNk7WL24iPWFcuz35ms//6fz5N+PAAApAJrF4NifXfv3776Z5KPBACA1GHtYlCs7+4tzr8uyUcCAEDqsHYxKJYbxgX9I/T90I3cWI6GhoaGxpVNWq9d6uvrdeeddyo3N1fjxo3T3Llz1dTUNKDp6elRdXW1xowZo9GjR2v+/Plqa2tL6EEnUiw3jKubPUU52VncWI6GhoaGxpVNWt9YrrGxUdXV1dq7d6+2b9+us2fP6uGHH1Z3d3e0qa2t1bZt27RlyxY1NjaqpaVF8+bNS/iBAwCAzDSktcvf//53jRs3To2NjbrvvvvU0dGhsWPHauPGjVqwYIEk6dixY5o8ebL27NmjGTNmXPU1Ta9dYrlhXCxT23S8URANDQ0NDY3rbix34sQJ3XrrrTpy5IimTp2qHTt2aObMmfr666+Vl5cX7YqLi1VTU6Pa2tpLXiMSiSgSiQw4+KKiIqa2AAAYlBFT276+PtXU1Ojuu+/W1KlTJUnhcFg5OTkDLjwkqaCgQOFw+LKvU19fL7/fH/0oKiq61kO6JkxtAQDIkKltdXW1jh49qs2bNw/pAFauXKmOjo7oR3Nz85BeL15MbQEAyICp7dKlS/Xee+9p586dGj9+fPTxQCCg3t5etbe3D+jb2toUCAQu+1per1c+n2/Ah0lMbWloaGhobG/SemrrOI6WLl2qrVu3aseOHSopKRnwfHl5uYYPH66GhoboY01NTTp16pRCoVBijjjBrnaDOUfnb7jD1JaGhoaGxq2N6altdjxxdXW1Nm7cqHfffVe5ubnR93H4/X6NHDlSfr9fixcv1vLly5Wfny+fz6dly5YpFArFtHRJlUemBvXEfSV6bffJyz7/+u6T+l8Trjd8VAAAuFNcaxeP5/JXRW+88YZ+8IMfSDr/S8ZWrFihTZs2KRKJqLKyUuvXr7/ij10uZnpqK8U2t023WRQNDQ0NDY2VU9tkSMXFR6xzWwAA3CojprZuYnJeBABAOsqIqa2bmJwXAQCQjkx+LYzrDadu1T+3DXf0XPbdwBf+zKytk4aGhoaGxl1NIJ2ntm4Vy51t020WRUNDQ0NDk6lTWy4+AACAUaxdxNSWhoaGhsbuhqktU1sAAIxjamsYU1sAgO2Y2hrG1BYAYDumtoYxtaWhoaGhsblhapsCTG1paGhoaGxumNoCAABXY+0iprY0NDQ0NHY3TG2Z2gIAYBxTW8OY2gIAbMfU1jCmtgAA2zG1NYypLQ0NDQ2NzQ1T2xRgaktDQ0NDY3PD1BYAALgaaxcxtaWhoaGhsbthasvUFgAA45jaGsbUFgBgO6a2hjG1BQDYjqmtYUxtaWhoaGhsbpjapgBTWxoaGhoamxumtgAAwNVYu4ipLQ0NDQ2N3Q1TW6a2AAAYx9TWMKa2AADbMbU1jKktAMB2Jr8WcvGhf09tBxPweRXw0dDQ0NDQuK8JMrU1b1iWR3NKg4M2VWWFqiqjoaGhoaFxXzOnNGh0assbTnX1tYt0/spxsHcL09DQ0NDQZGoTZO3C2gUAANNYuxjG2gUAYDvWLoaxdgEA2I4byxnGjeVoaGhoaGxuuLFcCnBjORoaGhoamxtuLAcAAFyNtYu4sRwNDQ0Njd0NN5ZjagsAgHFMbQ1jagsAsB1TW8OY2gIAbMfU1jCmtjQ0NDQ0NjdMbVOAqS0NDQ0Njc0NU1sAAOBqrF3E1JaGhoaGxu6GqS1TWwAAjGNqaxhTWwCA7ZjaGsbUFgBgO6a2hjG1paGhoaGxuWFqmwJMbWloaGhobG6Y2gIAAFdj7SKmtjQ0NDQ0djdMbZnaAgBgHFNbw5jaAgBsx9TWMKa2AADbMbU1jKktDQ0NDY3NTdpPbXfv3q3Zs2ersLBQHo9H77zzzoDnHcfRc889p2AwqJEjR2rWrFk6fvx4oo43KZja0tDQ0NDY3KT91La7u1ulpaVat27dZZ9/8cUXtXbtWr366qvat2+fRo0apcrKSvX08L4KAAAwxLWLx+PR1q1bNXfuXEnnv+tRWFioFStW6Omnn5YkdXR0qKCgQBs2bNDjjz9+1ddkaktDQ0NDQ8PU9sqffNHFx1//+lfdfPPNOnTokMrKyqLd/fffr7KyMr300kuXvEYkElEkEhlw8EVFRUxtAQAwKGOntuFwWJJUUFAw4PGCgoLocxerr6+X3++PfhQVFSXykGLC1BYAYDurprYrV65UR0dH9KO5udn4MTC1BQDYLmOntoFAQJLU1tamYDAYfbytrW3Aj2Eu5PV65fV6E3kYcWNqS0NDQ0Njc5P2U9vBlJSUKBAIqKGhIfpYZ2en9u3bp1AolMi/KqGY2tLQ0NDQ2NyYntrG/Z2PM2fO6MSJE9E/nzx5UocPH1Z+fr4mTJigmpoavfDCC7r11ltVUlKiZ599VoWFhdE3pQIAALvFvXbZtWuXHnjggUseX7RokTZs2CDHcVRXV6fXX39d7e3tuueee7R+/XrddtttMb0+U1saGhoaGhqmtkZxV1sAAMzL2KltpmJqCwCwnVVT23TA1BYAYLuMndpmKqa2NDQ0NDQ2Nxk9tc1UTG1paGhoaGxu0v6utgAAAEPB2kVMbWloaGho7G6Y2jK1BQDAOKa2hjG1BQDYjqmtYUxtAQC2M/m1kIsP/XtqO5iAz6uAj4aGhoaGxn1NkKmtecOyPJpTGhy0qSorVFUZDQ0NDQ2N+5o5pUGjU1vecKqrr12k81eOg71bmIaGhoaGJlObIGsX1i4AAJjG2sUw1i4AANuxdjGMtQsAwHbcWM4wbixHQ0NDQ2Nzw43lUoAby9HQ0NDQ2NxwYzkAAOBqrF3EjeVoaGhoaOxuuLEcU1sAAIxjamsYU1sAgO2Y2hrG1BYAYDumtoYxtaWhoaGhsblhapsCTG1paGhoaGxumNoCAABXY+0iprY0NDQ0NHY3TG2Z2gIAYBxTW8OY2gIAbMfU1jCmtgAA2zG1NYypLQ0NDQ2NzQ1T2xRgaktDQ0NDY3PD1BYAALgaaxcxtaWhoaGhsbthasvUFgAA45jaGsbUFgBgO6a2hjG1BQDYjqmtYUxtaWhoaGhsbpjapgBTWxoaGhoamxumtgAAwNVYu4ipLQ0NDQ2N3Q1TW6a2AAAYx9TWMKa2AADbMbU1jKktAMB2TG0NY2pLQ0NDQ2Nzw9Q2BZja0tDQ0NDY3DC1BQAArsbaRUxtaWhoaGjsbpjaMrUFAMA4praGMbUFANiOqa1hTG0BALZjamsYU1saGhoaGpsbprYpwNSWhoaGhsbmhqktAABwNdYuYmpLQ0NDQ2N3w9SWqS0AAMYxtTWMqS0AwHZMbQ1jagsAsJ3Jr4VJu/hYt26dbrzxRo0YMUIVFRX66KOPkvVXDVn/1HYwAZ9XAR8NDQ0NDY37mqAbprZvvvmmli9frrq6Oh08eFClpaWqrKzU6dOnk/HXDdmwLI/mlAYHbarKClVVRkNDQ0ND475mTmnQ6NQ2KW84raio0J133qmXX35ZktTX16eioiItW7ZMzzzzzKCfm45rF+n8leNg7xamoaGhoaHJ1CaY6WuX3t5eXXfddXr77bc1d+7c6OOLFi1Se3u73n333QF9JBJRJBIZcPBFRUWsXQAAMCij1y5ffvmlzp07p4KCggGPFxQUKBwOX9LX19fL7/dHP4qKihJ9SFfF2gUAYDur1i4rV65UR0dH9KO5udn4MbB2AQDYLqNvLHfDDTdo2LBhamtrG/B4W1ubAoHAJb3X65XX6030YcSFG8vR0NDQ0NjcZPyN5XJyclReXq6GhoboY319fWpoaFAoFEr0X5cQ/TeWk87/R7hQ/5/7b8xDQ0NDQ0PjtsYVN5Zbvny5fvWrX+m3v/2tPv/8cy1ZskTd3d364Q9/mIy/LiEemRrUK9+brsBFv+8j4B+hV743XY9MDdLQ0NDQ0Li2MSlp93Z5+eWXtWbNGoXDYZWVlWnt2rWqqKi46uelYmp7oXN9jj46+ZVOd/VoXO75b0NdfDVIQ0NDQ0Pj1uZacWM5AABgFDeWAwAAaYuLDwAAYBQXHwAAwCguPgAAgFFcfAAAAKO4+AAAAEZx8QEAAIzi4gMAABjFxQcAADAq4Xe1Har+X7ja2dmZ4iMBAACx6v+6HcsvTk+7i4+uri5JUlFRUYqPBAAAxKurq0t+v3/QJu3u7dLX16eWlhbl5ubK4zF3e98LdXZ2qqioSM3NzdxfJok4z2Zwns3gPJvBeTbjWs6z4zjq6upSYWGhsrIGf1dH2n3nIysrS+PHj0/1YUiSfD4f/7gN4DybwXk2g/NsBufZjHjP89W+49GPN5wCAACjuPgAAABGcfFxGV6vV3V1dfJ6vak+FFfjPJvBeTaD82wG59mMZJ/ntHvDKQAAcDe+8wEAAIzi4gMAABjFxQcAADCKiw8AAGAUFx8XWbdunW688UaNGDFCFRUV+uijj1J9SBlt9+7dmj17tgoLC+XxePTOO+8MeN5xHD333HMKBoMaOXKkZs2apePHj6fmYDNYfX297rzzTuXm5mrcuHGaO3eumpqaBjQ9PT2qrq7WmDFjNHr0aM2fP19tbW0pOuLM9Morr2jatGnRX7wUCoX0/vvvR5/nHCfH6tWr5fF4VFNTE32Mcz10zz//vDwez4CPSZMmRZ9P5jnm4uMCb775ppYvX666ujodPHhQpaWlqqys1OnTp1N9aBmru7tbpaWlWrdu3WWff/HFF7V27Vq9+uqr2rdvn0aNGqXKykr19PQYPtLM1tjYqOrqau3du1fbt2/X2bNn9fDDD6u7uzva1NbWatu2bdqyZYsaGxvV0tKiefPmpfCoM8/48eO1evVqHThwQB9//LEefPBBVVVV6dNPP5XEOU6G/fv367XXXtO0adMGPM65Tozbb79dra2t0Y8PP/ww+lxSz7GDqLvuusuprq6O/vncuXNOYWGhU19fn8Kjcg9JztatW6N/7uvrcwKBgLNmzZroY+3t7Y7X63U2bdqUgiN0j9OnTzuSnMbGRsdxzp/X4cOHO1u2bIk2n3/+uSPJ2bNnT6oO0xWuv/5659e//jXnOAm6urqcW2+91dm+fbtz//33O0899ZTjOPx7TpS6ujqntLT0ss8l+xzznY9/6e3t1YEDBzRr1qzoY1lZWZo1a5b27NmTwiNzr5MnTyocDg84536/XxUVFZzzIero6JAk5efnS5IOHDigs2fPDjjXkyZN0oQJEzjX1+jcuXPavHmzuru7FQqFOMdJUF1dre985zsDzqnEv+dEOn78uAoLC3XTTTdp4cKFOnXqlKTkn+O0u7Fcqnz55Zc6d+6cCgoKBjxeUFCgY8eOpeio3C0cDkvSZc95/3OIX19fn2pqanT33Xdr6tSpks6f65ycHOXl5Q1oOdfxO3LkiEKhkHp6ejR69Ght3bpVU6ZM0eHDhznHCbR582YdPHhQ+/fvv+Q5/j0nRkVFhTZs2KCJEyeqtbVVq1at0r333qujR48m/Rxz8QG4THV1tY4ePTrgZ7dInIkTJ+rw4cPq6OjQ22+/rUWLFqmxsTHVh+Uqzc3Neuqpp7R9+3aNGDEi1YfjWo8++mj0f0+bNk0VFRUqLi7WW2+9pZEjRyb17+bHLv9yww03aNiwYZe8k7etrU2BQCBFR+Vu/eeVc544S5cu1XvvvaedO3dq/Pjx0ccDgYB6e3vV3t4+oOdcxy8nJ0e33HKLysvLVV9fr9LSUr300kuc4wQ6cOCATp8+renTpys7O1vZ2dlqbGzU2rVrlZ2drYKCAs51EuTl5em2227TiRMnkv7vmYuPf8nJyVF5ebkaGhqij/X19amhoUGhUCiFR+ZeJSUlCgQCA855Z2en9u3bxzmPk+M4Wrp0qbZu3aodO3aopKRkwPPl5eUaPnz4gHPd1NSkU6dOca6HqK+vT5FIhHOcQDNnztSRI0d0+PDh6Mcdd9yhhQsXRv835zrxzpw5oy+++ELBYDD5/56H/JZVF9m8ebPj9XqdDRs2OJ999pnzxBNPOHl5eU44HE71oWWsrq4u59ChQ86hQ4ccSc7PfvYz59ChQ87f/vY3x3EcZ/Xq1U5eXp7z7rvvOp988olTVVXllJSUON98802KjzyzLFmyxPH7/c6uXbuc1tbW6Mc///nPaPPkk086EyZMcHbs2OF8/PHHTigUckKhUAqPOvM888wzTmNjo3Py5Ennk08+cZ555hnH4/E4f/rTnxzH4Rwn04VrF8fhXCfCihUrnF27djknT550/vznPzuzZs1ybrjhBuf06dOO4yT3HHPxcZFf/vKXzoQJE5ycnBznrrvucvbu3ZvqQ8poO3fudCRd8rFo0SLHcc7PbZ999lmnoKDA8Xq9zsyZM52mpqbUHnQGutw5luS88cYb0eabb75xfvzjHzvXX3+9c9111znf/e53ndbW1tQddAb60Y9+5BQXFzs5OTnO2LFjnZkzZ0YvPByHc5xMF198cK6H7rHHHnOCwaCTk5PjfOtb33Iee+wx58SJE9Hnk3mOPY7jOEP//gkAAEBseM8HAAAwiosPAABgFBcfAADAKC4+AACAUVx8AAAAo7j4AAAARnHxAQAAjOLiAwAAGMXFBwAAMIqLDwAAYBQXHwAAwCguPgAAgFH/H451iHU40DlbAAAAAElFTkSuQmCC",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def neighbors(elm: ELM, row: int, col: int) -> Iterator[Tuple[float, int, int]]:\n",
|
|
" for row_offset in [-1, 0, 1]:\n",
|
|
" for col_offset in [-1, 0, 1]:\n",
|
|
" nb_row = row + row_offset\n",
|
|
" nb_col = col + col_offset\n",
|
|
" \n",
|
|
" if nb_row < 0 or nb_row >= elm.n_rows or nb_col < 0 or nb_col >= elm.n_cols:\n",
|
|
" continue\n",
|
|
" if row_offset == 0 and col_offset == 0:\n",
|
|
" continue\n",
|
|
" \n",
|
|
" k = 1.0 if row_offset == 0 or col_offset == 0 else 2.0\n",
|
|
" yield (k, nb_row, nb_col)\n",
|
|
"\n",
|
|
"def calculate_force(elm: ELM, displacement: numpy.ndarray, row: int, col: int) -> numpy.ndarray:\n",
|
|
" K = elm.elastic_spring_constant\n",
|
|
" c = elm.bond_bending_constant\n",
|
|
"\n",
|
|
" total = numpy.array([0.0, 0.0])\n",
|
|
" for k, nb_row, nb_col in neighbors(elm, row, col):\n",
|
|
" u_ij = displacement[row, col, :] - displacement[nb_row, nb_col, :]\n",
|
|
" x_ij = elm.get_lattice_position(row, col) - elm.get_lattice_position(nb_row, nb_col)\n",
|
|
" r_ij = x_ij + u_ij\n",
|
|
" \n",
|
|
" total -= K * k * numpy.dot(u_ij, x_ij) * (r_ij / numpy.linalg.norm(r_ij)) + (c * u_ij) / (numpy.linalg.norm(x_ij) ** 2)\n",
|
|
"\n",
|
|
" return total\n",
|
|
"\n",
|
|
"def gaussian(mu, sigma_sq, x):\n",
|
|
" return numpy.exp(-(((x - mu) / numpy.sqrt(sigma_sq)) ** 2) / 2) / numpy.sqrt(2 * numpy.pi * sigma_sq)\n",
|
|
"\n",
|
|
"e = ELM.make(50, 50)\n",
|
|
"n_plucked = 10\n",
|
|
"base = e.n_rows // 2 - n_plucked // 2\n",
|
|
"for row in range(n_plucked):\n",
|
|
" e.set_displacement(base + row, 0, numpy.array([\n",
|
|
" gaussian(n_plucked // 2, n_plucked / 3, row) * -5,\n",
|
|
" 0.0\n",
|
|
" ]))\n",
|
|
"\n",
|
|
"x, y = numpy.meshgrid(\n",
|
|
" numpy.linspace(0, e.n_rows - 1, e.n_rows, dtype=int),\n",
|
|
" numpy.linspace(0, e.n_cols - 1, e.n_cols, dtype=int)\n",
|
|
")\n",
|
|
"\n",
|
|
"plt.scatter(\n",
|
|
" x + e.displacement[:, :, 0],\n",
|
|
" y + e.displacement[:, :, 1]\n",
|
|
")\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "8dc183cc-5fa4-48f1-bc30-9cffff91be6a",
|
|
"metadata": {},
|
|
"source": [
|
|
"A [velocity Verlet integration](https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet) scheme is used to solve the system.\n",
|
|
"\n",
|
|
"Our initial configuration defines for each of the $M$ nodes $i$:\n",
|
|
"\n",
|
|
"* $\\mathbf{u}_i(t_0) \\in \\mathbb{R}^M \\times 2$: Displacement\n",
|
|
"* $\\mathbf{v}_i(t_0) \\in \\mathbb{R}^M \\times 2$: Velocity\n",
|
|
"\n",
|
|
"To calculate a time step, do for every node $i$:\n",
|
|
"\n",
|
|
"* Calculate the displacement $\\mathbf{u}_i(t + \\Delta t) = \\mathbf{u}_i(t) + \\mathbf{v}_i(t) \\Delta t + \\frac{1}{2} \\mathbf{F}_i(t) \\Delta t^2$\n",
|
|
"* Calculate $\\mathbf{F}_i(t + \\Delta t)$ using $\\mathbf{u}_i(t + \\Delta t)$\n",
|
|
"* Calculate $\\mathbf{v}_i(t + \\Delta t) = \\mathbf{v}_i(t) + \\frac{1}{2} (\\mathbf{F}_i(t) + \\mathbf{F}_i(t + \\Delta t)) \\Delta t$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "ff7adc30-191e-4b7c-b146-2a533463d6ea",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def velocity_verlet_step(elm: ELM, delta_t: float) -> ELM:\n",
|
|
" new_displacement = elm.displacement + elm.velocity * delta_t + 0.5 * elm.force * delta_t * delta_t\n",
|
|
" \n",
|
|
" new_force = numpy.zeros((elm.n_rows, elm.n_cols, 2))\n",
|
|
" for row in range(elm.n_rows):\n",
|
|
" for col in range(elm.n_cols):\n",
|
|
" new_force[row, col, :] = calculate_force(\n",
|
|
" elm,\n",
|
|
" new_displacement,\n",
|
|
" row,\n",
|
|
" col\n",
|
|
" )\n",
|
|
" \n",
|
|
" new_velocity = elm.velocity + 0.5 * (elm.force + new_force) * delta_t\n",
|
|
"\n",
|
|
" return ELM(elm.n_rows, elm.n_cols, new_displacement, new_velocity, new_force)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"id": "bea964bf-6e94-4b72-a458-50c2964bbba1",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"............................................................................................................................................................................................................................................................................................................"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"es = [e]\n",
|
|
"for _ in range(300):\n",
|
|
" print(\".\", end=\"\")\n",
|
|
" es.append(velocity_verlet_step(es[-1], 0.1))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"id": "421c0fd2-2dba-4fd5-89ff-0968227b08b9",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "dc96f81e9f764b4fa58103f705f3e445",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"interactive(children=(IntSlider(value=0, description='i', max=300), Output()), _dom_classes=('widget-interact'…"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<function __main__.do_plot(i: int) -> None>"
|
|
]
|
|
},
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"def do_plot(i: int) -> None:\n",
|
|
" fig, ax = plt.subplots(figsize=(10, 5), ncols=2)\n",
|
|
" ax[0].scatter(\n",
|
|
" x + es[i].displacement[:, :, 0],\n",
|
|
" y + es[i].displacement[:, :, 1]\n",
|
|
" )\n",
|
|
" ax[0].quiver(\n",
|
|
" x + es[i].displacement[:, :, 0],\n",
|
|
" y + es[i].displacement[:, :, 1],\n",
|
|
" es[i].force[:, :, 0],\n",
|
|
" es[i].force[:, :, 1],\n",
|
|
" angles = \"xy\",\n",
|
|
" scale = 100\n",
|
|
" )\n",
|
|
" ax[1].imshow(numpy.linalg.norm(es[i].displacement, axis=2), cmap=\"hot\")\n",
|
|
" plt.show()\n",
|
|
"\n",
|
|
"interact(do_plot, i = IntSlider(min=0, max=len(es) - 1, value=0))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 41,
|
|
"id": "96ab8cf1-8739-4faa-b6b5-70aa52a15a2a",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAG4CAYAAAAE8e26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAl3klEQVR4nO3df3RU5Z3H8U9CyCSQzEAC5kdJkFPQKBSsQeJody0QfoirIqHV1m5Zyx4PGqmYUmy6i7QVTlA4KrgoWDxo9xRx4xEtnC0sjRp1DRTiRhAki126pA0J0G4yIZJJTJ79gzJ1SNDMZPLcSfJ+nTNHc3/Nd577ZT5z78ydiTHGGAEAAGtinS4AAICBhvAFAMAywhcAAMsIXwAALCN8AQCwjPAFAMAywhcAAMsIXwAALCN8AQCwjPAFAMCyOKcL6Kt8Pp9qa2vDWnfcuHEaNGhQhCtCtKNnECp6pv8ifMP0q1/9Sn//938f1rqnT5/WiBEjIlwRoh09g1DRM/0Xp53DNHjwYKdLQB9DzyBU9Ez/RfiGqaCgQLfeemvg7x/96EcyxnTrxqvRgYmeQajomX7MhKGpqcnExMSYJ554IpzV+41z586Zr3/960aSkWQee+yxkLcxevRos2LFiojUE8ltoXfQMwiVEz3T0dFhhg4dapYtWxaR7aGzsI58P/zwQxljNHHixMi9CriI3+/Xww8/rMzMTCUmJiovL0979uzptfsLR0JCgnbs2KEpU6ZIkh5++GFt2rQpItvevXu3YmJiurzdeeedEbmPcPV034SyfneX7Qv9IvVuz0QzeiZ8TvTM73//ezU3N+srX/lKr97P5zl79qxWrFih2bNnKyUlRTExMXrhhRdC2kZlZaVmz54tt9ut5ORkzZw5U1VVVV0u+/777+u2225TSkqKhgwZogkTJmj9+vVBy7z11luXfF7eu3dvaA8wnMRub283586dMx0dHRF+LfBXd911l4mLizNLly41mzZtMl6v18TFxZl33nmn1+4zXH/605/MhAkTjCQTGxtrXnrppW6ve6lXkI899piRZNavX2/+9V//NehWVVUV0rYiraf7JpT1u7tsX+oXY3qnZ8JBz9AzXdmxY4eR5OhzzfHjx40kk52dHTjy37JlS7fXr6ysNAkJCWbcuHFm7dq15vHHHzeXX365cbvd5ujRo0HL7t6928THx5u8vDzzxBNPmOeee848/PDD5oc//GHQcm+++aaRZL7//e93el4+ffp0SI8vrPDtbfv27TOSzJo1awLTzp07Z7785S8br9frYGWXVltba7785S8bSWbw4MFm586d3VrvUk189913G4/HE9ILHBv/IHq6b0JZv7vL9sV+MSbyPRMOeoae6crq1atNXFyc8fv9EdleOFpaWszJkyeNMcbs378/5PCdM2eOGT58uDlz5kxgWm1trUlKSjLz5s0LTGtsbDRpaWnmjjvuMO3t7Z+7zQvhW1paGtqD6UJYp51nzJihG2+8MZxVu+WVV17RoEGDdO+99wamJSQkaOHChaqoqFBNTU2v3Xe4MjIy9Jvf/EZf+tKX1NbWpvnz56u8vDzs7X3wwQf66le/qpiYmAhW2XM93TehrN/dZftiv0iR7ZkZM2Zo7NixXc6bPHmy8vLyelJqj9AzkRPp5xlJevnll3XNNdcoISFBubm5+u1vf6vDhw/riiuuUHx8fIQqD53L5VJ6enrY67/zzjvKz89XampqYFpGRoZuuukm7dy5U2fPnpUkbd26VfX19Vq1apViY2PV3Nysjo6OL9x+U1OTPv3007DrCyt8Dx482OX7vW1tbTpz5ky3bp/34P7rv/5LV1xxhdxud9D0C+95XOqcvdMuv/xy7dmzRyNGjFBLS4tuvfVWHThwIOTttLa2qrq6Wjk5OZ3Gra2tLez6IrF/erpvQlm/u8v21X6RItcz48eP1+9//3v5/f6g6a+99poqKyu1cuXKsOqjZ6JPpHpGkp588kndddddGjNmjNavXy+v16u/+7u/0969e3v0fm+ksqAn/H6/EhMTO00fMmSIWltb9eGHH0qSfvOb38jtduuPf/yjrrzySiUlJcntduu+++5TS0tLl9u+55575Ha7lZCQoKlTp4Y1/iF/ycapU6d06tSpLsP3P//zPzV16tRubef48eO6/PLLu5x38uRJZWRkdJp+YVq43/hiw1VXXaXdu3dr6tSp8vl8mj17tt5++21dffXV3d7GkSNH1NbWpo0bN2rjxo1B8w4fPhzStj4rEvunp/smlPW7u2xf7hcpMj0zfvx4tbe369ixY5owYYIkyRijFStW6KabbtKMGTPCqo2eiU6R6JmqqiotW7ZMP/7xj7Vq1arA9I6ODj377LP67ne/G3Z9kcqCnrjyyiu1d+9etbe3B77pq7W1Vfv27ZMk/fGPf5QkHTt2TJ9++qluv/12LVy4UCUlJXrrrbf09NNPq6GhQS+99FJgm/Hx8SooKNCcOXM0YsQIHTlyRGvXrtXf/M3f6L333tNXv/rVbtcXcvgePHhQkrp8VTRp0qRuf1rw804nnDt3Ti6Xq9P0hISEwPxodu2112rOnDnatm2b/vSnP+nVV18N6R/FhTF+8cUXlZmZGTQvJydH0vlXddnZ2Tp27FinV++XEon909N9E8r63V22r/eL1POeGT9+vCTp6NGjgfAtLS3VwYMH9c4770iiZz67LD0jrVq1Sh6PR//0T/8UNP2mm27Ss88+G3iOd6pveur+++/Xfffdp4ULF2rZsmXq6OjQypUrdfLkSUl/3cdnz57VJ598okWLFgU+3Txv3jy1trZq06ZN+tnPfqZx48ZJkm644QbdcMMNgfu47bbbNH/+fE2cOFHFxcXatWtXt+sLOXwPHTokqevwHT58uPLz80PdZCeJiYmdTp9JCpwC6OpUQjRZuXKltm3bJklasGCB/vmf/zmk9T/44APFx8fr29/+tuLiut5FLpdL9fX1IW03Evunp/smlPW7u2xf7xep5z3z2fCVzh+9/OQnP9GsWbP0ta99TRI989llB3rP+P1+/fu//7vuvfdeDRkyJGjehfcxLzzHO9U3PbVo0SLV1NRozZo1evHFFyWd//zDsmXLtGrVKiUlJUn6677+1re+FbT+t7/9bW3atEkVFRWB8O3K2LFjdfvtt+vVV18NOsr+ImEd+WZnZ8vj8XSa19raqj//+c/d2s7IkSMvWWRGRkbglMBnXXjFcvHRYDTZvHmzli9fLkm65ZZbtHnz5pC3cfDgQY0bN+6SwRuuSOyfnu6bUNbv7rJ9uV+kyPSMx+NRZmZmIHy3bt2qjz76SL/4xS96VBs9E5162jO/+93v9Mknnyg3N7fTvAMHDigpKUljxowJu75IZUFPrVq1SkuXLtXhw4fl8Xj0la98RT/+8Y8lSVdccYWk8/v68OHDSktLC1r3sssukyT93//93xfeT1ZWllpbW9Xc3NztswNhhe+lvlzjvffei8h5/muuuUZvvvmmfD5f0AO5cK7+mmuuCalmW15//XUtWrRI0vnTE6WlpWEF6MGDB/X1r3/9c5dZt26dDh48qOeff77b243E/unpvgll/e4u21f7RYpcz0jnj36PHj2q9vZ2/exnP9PcuXM1efLkwHx6hp654JNPPulyenNzs37xi19o/PjxgSstnOqbSBk+fHjg7I90/gNWo0aNCryFl5ubqz179gQ+cHXBhff9R44c+YX38T//8z9KSEgIHE13R0h7rL29XUeOHNHs2bO7nB+p8/zz58/X2rVr9dxzz2np0qWSzp8m2bJli/Ly8pSVlRVK2Va8++67uuuuu9Te3q7x48dr586dYZ26qqur06lTpwKnES/l814EXUok9k8o++aTTz7RiRMnNGLEiMD3zIayfneX7Yv9IkWuZy4YP368Nm/erBdffFG/+93v9OqrrwbNp2fomQtGjx4tSXrjjTf0ne98JzB95cqV+vOf/xz0tqJTfROKrvqmKy+//LL279+vtWvXKjb2/MU+3/zmN7V69Wo9//zzmjZtWmDZzZs3Ky4uLuhA6PTp053C+IMPPtCvfvUr3XzzzYFtdksoFwV/9NFHRlJI36wSrm984xsmLi7O/PCHPzSbNm0yN9xwg4mLizPl5eW9ft+hOnTokBk2bJiRZLKyskxNTU231734YvVdu3YZSeaVV1753PUmT55s3njjjc/dVm/p7r65cEH6xTWFsm+7u2xf6hdjItszF/z85z83kszIkSPNt771rU7z6Rl65rNmzpxpYmJizKJFi8ymTZvMHXfcYUaOHGkkmXXr1gWWc7Jvnn76afPoo4+a++67z0gy8+bNM48++qh59NFHTUNDQ2C5rvqmvLzcTJ8+3Tz22GNm8+bN5h//8R/NoEGDzOzZs01bW1vQ/Xzve98zksw3v/lNs2HDBvONb3zDSDLFxcVBy02dOtXMmTPHrFy50jz33HNmyZIlZsiQIcbj8ZgjR46E9NhCCt9/+7d/M5LM4cOHQ7qTcJw7d84sXbrUpKenG5fLZa677jqza9euXr/fUP3v//6v+dKXvmQkmdTUVPPRRx+FtP7FTfz4448bSZ+7I9vb282QIUOCvrmlq231lu7um0s9kYayb7u7bF/pF2Mi3zMXvPfee0aSGTRokPnv//7voHn0DD1z8ZiePHnS3HbbbSY5OdmkpqaaO++80/zyl780kkxZWZkxxvm+GT16dOAHJS6+HT9+PLBcV33z8ccfm5kzZ5oRI0YYl8tlcnJyTElJSZff2tXa2mp+8pOfmNGjR5vBgwebsWPHmieffLLTcuvWrTNTpkwxKSkpJi4uzmRkZJjvfOc75tixYyE/tqj8esm+4syZMyYnJ8dIMkOHDjV79+4NeRvhNHF1dbXJyMiIyLZgFz2DUDnVM8bQN72J3/MNU3Nzs2655RYdPXpUcXFxKi0ttfYVfuG8BwPn0TMIlZM9I9E3vYnwDdODDz4Y+GTkp59+qjlz5lzyp6YuvoX6s1gXO3ToEP8g+iB6BqFysmck+qY3RfZC0gHkD3/4g2P3/dOf/tSx+0b46BmEysmekeib3hRjjDFOF9EX/eEPfwj8KkaoMjIyuvySEvRv9AxCRc/0X4QvAACW8Z4vAACW9dv3fDs6OlRbW6vk5OSo+0F6AMAXM8aoqalJmZmZoX17VB8Q1eG7YcMGrVmzRnV1dZo0aZKefvrpwI9df5Ha2tqo/Yo4AED31dTUaNSoUU6XEVFRG74vv/yyioqKtHHjRuXl5empp57SrFmzVF1dHfi1ic+TnJwsSUqQxHEvAPQ9RlKL/vp83p9E7Qeu8vLydN111+lf/uVfJJ0/jZyVlaXFixfrRz/6Uafl/X5/0O9z+nw+ZWVlKVGELwD0RUbSOUmNjY3d/qm+viIqT6K3traqsrIy6MeYY2NjlZ+fr4qKii7XKSkpkcfjCdw45QwAiFZRGb5nzpxRe3t7px83TktLU11dXZfrFBcXq7GxMXCrqamxUSoAACGL2vd8Q+VyueRyuZwuAwCALxSVR74jRozQoEGDVF9fHzS9vr4+Ij+8DACAk6IyfOPj45Wbm6uysrLAtI6ODpWVlcnr9TpYGQAAPRe1p52Lioq0YMECTZ48WVOmTNFTTz2l5uZm3XPPPU6XBgBAj0Rt+N555506ffq0HnnkEdXV1emaa67Rrl27On0ICwCAviZqr/PtKZ/PJ4/Hw3W+ANBHcZ0vAACIGMIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLCF8AACxzJHzffvtt3XrrrcrMzFRMTIxee+21oPnGGD3yyCPKyMhQYmKi8vPzdezYMSdKBQAg4hwJ3+bmZk2aNEkbNmzocv7jjz+u9evXa+PGjdq3b5+GDh2qWbNmqaWlxXKlAABEXowxxjhaQEyMtm/frrlz50o6f9SbmZmpH/zgB1q6dKkkqbGxUWlpaXrhhRd01113dbkdv98vv98f+Nvn8ykrK0uJkmJ6+0EAACLOSDqn8xngdrudLieiou493+PHj6uurk75+fmBaR6PR3l5eaqoqLjkeiUlJfJ4PIFbVlaWjXIBAAhZ1IVvXV2dJCktLS1oelpaWmBeV4qLi9XY2Bi41dTU9GqdAACEK87pAiLF5XLJ5XI5XQYAAF8o6o5809PTJUn19fVB0+vr6wPzAADoy6IufMeMGaP09HSVlZUFpvl8Pu3bt09er9fBygAAiAxHTjufPXtWH3/8ceDv48ePq6qqSikpKcrOztaSJUu0cuVKjRs3TmPGjNHy5cuVmZkZ+EQ0AAB9mSPhe+DAAU2dOjXwd1FRkSRpwYIFeuGFF7Rs2TI1Nzfr3nvvVUNDg772ta9p165dSkhIcKJcAAAiyvHrfHuLz+eTx+PhOl8A6KO4zhcAAEQM4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGWELwAAlhG+AABYRvgCAGAZ4QsAgGVxThfQ2zzq+hVGg+U6AAC4wJEj35KSEl133XVKTk7WZZddprlz56q6ujpomZaWFhUWFio1NVVJSUkqKChQfX29E+UCABBRjoRveXm5CgsLtXfvXu3Zs0dtbW2aOXOmmpubA8s89NBD2rFjh0pLS1VeXq7a2lrNmzfPiXIBAIioGGOMcbqI06dP67LLLlN5ebn+9m//Vo2NjRo5cqS2bt2q+fPnS5KOHj2qq666ShUVFbr++us7bcPv98vv9wf+9vl8ysrKUro47QwAfZGRdE5SY2Oj3G630+VEVFR84KqxsVGSlJKSIkmqrKxUW1ub8vPzA8vk5OQoOztbFRUVXW6jpKREHo8ncMvKyur9wgEACIPj4dvR0aElS5boxhtv1IQJEyRJdXV1io+P17Bhw4KWTUtLU11dXZfbKS4uVmNjY+BWU1PT26UDABAWxz/tXFhYqA8//FDvvvtuj7bjcrnkcrkiVBUAAL3H0fB94IEHtHPnTr399tsaNWpUYHp6erpaW1vV0NAQdPRbX1+v9PT0kO7janX9IHsW9QAAhM+R087GGD3wwAPavn273njjDY0ZMyZofm5urgYPHqyysrLAtOrqap04cUJer9d2uQAARJQjR76FhYXaunWrXn/9dSUnJwfex/V4PEpMTJTH49HChQtVVFSklJQUud1uLV68WF6vt8tPOgMA0Jc4cqlRTExMl9O3bNmif/iHf5B0/ks2fvCDH+ill16S3+/XrFmz9Mwzz3T7tLPP55PH49E0cdoZAPqi/nypUVRc59sbCF8A6Nv6c/g6fqkRAAADDeELAIBlhC8AAJY5/iUbve1+SUO7mM57vgAAp3DkCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACW9ftLjWYsldxd/czvKuulAAAgiSNfAACsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwLIYY4xxuoje4PP55PF41Nj4fbm7uNB3aMwaB6oCAHSXkXROUmNjo9xut9PlRBRHvgAAWEb4AgBgGeELAIBlhC8AAJYRvgAAWEb4AgBgWb//SUHpUUldfUSdS40AAM7gyBcAAMsIXwAALCN8AQCwjPAFAMAywhcAAMsIXwAALBsAlxotl9T5V40AAHAKR74AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZf3/Ot9H13OZLwAgqnDkCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACW9ftLjfaslYY6XQQAAJ/BkS8AAJYRvgAAWEb4AgBgGeELAIBlhC8AAJYRvgAAWEb4AgBgWb+/zvcZDYAHCQDoUzjyBQDAMsIXAADLCF8AACwjfAEAsIzwBQDAMsIXAADLHAnfZ599VhMnTpTb7Zbb7ZbX69Wvf/3rwPyWlhYVFhYqNTVVSUlJKigoUH19fVj3VSHp3S5uAAA4xZHwHTVqlFavXq3KykodOHBA06ZN0+23367Dhw9Lkh566CHt2LFDpaWlKi8vV21trebNm+dEqQAARFyMMcY4XYQkpaSkaM2aNZo/f75GjhyprVu3av78+ZKko0eP6qqrrlJFRYWuv/76bm3P5/PJ4/EoUVJML9YNAOgdRtI5SY2NjXK73U6XE1GOv+fb3t6ubdu2qbm5WV6vV5WVlWpra1N+fn5gmZycHGVnZ6uiouKS2/H7/fL5fEE3AACikWPhe+jQISUlJcnlcmnRokXavn27rr76atXV1Sk+Pl7Dhg0LWj4tLU11dXWX3F5JSYk8Hk/glpWV1cuPAACA8DgWvldeeaWqqqq0b98+3XfffVqwYIGOHDkS9vaKi4vV2NgYuNXU1ESwWgAAIsex3xyIj4/X2LFjJUm5ubnav3+/1q1bpzvvvFOtra1qaGgIOvqtr69Xenr6Jbfncrnkcrl6u2wAAHrM8fd8L+jo6JDf71dubq4GDx6ssrKywLzq6mqdOHFCXq/XwQoBAIgMR458i4uLdfPNNys7O1tNTU3aunWr3nrrLe3evVsej0cLFy5UUVGRUlJS5Ha7tXjxYnm93m5/0hkAgGjmSPieOnVK3/3ud3Xy5El5PB5NnDhRu3fv1owZMyRJTz75pGJjY1VQUCC/369Zs2bpmWeecaJUAAAiLmqu8400rvMFgL6N63wBAEDEEL4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYBnhCwCAZYQvAACWEb4AAFhG+AIAYJnj4bt69WrFxMRoyZIlgWktLS0qLCxUamqqkpKSVFBQoPr6eueKBAAgghwN3/3792vTpk2aOHFi0PSHHnpIO3bsUGlpqcrLy1VbW6t58+Y5VCUAAJHlWPiePXtWd999t37+859r+PDhgemNjY16/vnn9cQTT2jatGnKzc3Vli1b9N5772nv3r1OlQsAQMQ4Fr6FhYW65ZZblJ+fHzS9srJSbW1tQdNzcnKUnZ2tioqKS27P7/fL5/MF3QAAiEZxTtzptm3b9P7772v//v2d5tXV1Sk+Pl7Dhg0Lmp6Wlqa6urpLbrOkpEQ//elPI10qAAARZ/3It6amRg8++KB++ctfKiEhIWLbLS4uVmNjY+BWU1MTsW0DABBJ1sO3srJSp06d0rXXXqu4uDjFxcWpvLxc69evV1xcnNLS0tTa2qqGhoag9err65Wenn7J7bpcLrnd7qAbAADRyPpp5+nTp+vQoUNB0+655x7l5OTo4YcfVlZWlgYPHqyysjIVFBRIkqqrq3XixAl5vV7b5QIAEHHWwzc5OVkTJkwImjZ06FClpqYGpi9cuFBFRUVKSUmR2+3W4sWL5fV6df3119suFwCAiHPkA1df5Mknn1RsbKwKCgrk9/s1a9YsPfPMM06XBQBARMQYY4zTRfQGn88nj8ejREkxThcDAAiZkXRO57//ob99jsfxr5cEAGCgIXwBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALCM8AUAwDLCFwAAywhfAAAsI3wBALAszukCeosx5vx/Ha4DABCeC8/fF57P+5N+G75NTU2SpBaH6wAA9ExTU5M8Ho/TZURUjOmPLykkdXR0qLa2VsnJyYqJiZHP51NWVpZqamrkdrudLi9qMU7dwzh9McaoexinSzPGqKmpSZmZmYqN7V/vkvbbI9/Y2FiNGjWq03S3202DdwPj1D2M0xdjjLqHcepafzvivaB/vZQAAKAPIHwBALBswISvy+XSihUr5HK5nC4lqjFO3cM4fTHGqHsYp4Gp337gCgCAaDVgjnwBAIgWhC8AAJYRvgAAWEb4AgBgGeELAIBlAyJ8N2zYoMsvv1wJCQnKy8vTb3/7W6dLctTbb7+tW2+9VZmZmYqJidFrr70WNN8Yo0ceeUQZGRlKTExUfn6+jh075kyxDiopKdF1112n5ORkXXbZZZo7d66qq6uDlmlpaVFhYaFSU1OVlJSkgoIC1dfXO1SxM5599llNnDgx8A1NXq9Xv/71rwPzGaPOVq9erZiYGC1ZsiQwjXEaWPp9+L788ssqKirSihUr9P7772vSpEmaNWuWTp065XRpjmlubtakSZO0YcOGLuc//vjjWr9+vTZu3Kh9+/Zp6NChmjVrllpaBtbPVJSXl6uwsFB79+7Vnj171NbWppkzZ6q5uTmwzEMPPaQdO3aotLRU5eXlqq2t1bx58xys2r5Ro0Zp9erVqqys1IEDBzRt2jTdfvvtOnz4sCTG6GL79+/Xpk2bNHHixKDpjNMAY/q5KVOmmMLCwsDf7e3tJjMz05SUlDhYVfSQZLZv3x74u6Ojw6Snp5s1a9YEpjU0NBiXy2VeeuklByqMHqdOnTKSTHl5uTHm/LgMHjzYlJaWBpb56KOPjCRTUVHhVJlRYfjw4Wbz5s2M0UWamprMuHHjzJ49e8xNN91kHnzwQWMMvTQQ9esj39bWVlVWVio/Pz8wLTY2Vvn5+aqoqHCwsuh1/Phx1dXVBY2Zx+NRXl7egB+zxsZGSVJKSookqbKyUm1tbUFjlZOTo+zs7AE7Vu3t7dq2bZuam5vl9XoZo4sUFhbqlltuCRoPiV4aiPrtrxpJ0pkzZ9Te3q60tLSg6WlpaTp69KhDVUW3uro6SepyzC7MG4g6Ojq0ZMkS3XjjjZowYYKk82MVHx+vYcOGBS07EMfq0KFD8nq9amlpUVJSkrZv366rr75aVVVVjNFfbNu2Te+//77279/faR69NPD06/AFIqWwsFAffvih3n33XadLiUpXXnmlqqqq1NjYqFdeeUULFixQeXm502VFjZqaGj344IPas2ePEhISnC4HUaBfn3YeMWKEBg0a1OkTg/X19UpPT3eoquh2YVwYs7964IEHtHPnTr355ptBvxGdnp6u1tZWNTQ0BC0/EMcqPj5eY8eOVW5urkpKSjRp0iStW7eOMfqLyspKnTp1Stdee63i4uIUFxen8vJyrV+/XnFxcUpLS2OcBph+Hb7x8fHKzc1VWVlZYFpHR4fKysrk9XodrCx6jRkzRunp6UFj5vP5tG/fvgE3ZsYYPfDAA9q+fbveeOMNjRkzJmh+bm6uBg8eHDRW1dXVOnHixIAbq4t1dHTI7/czRn8xffp0HTp0SFVVVYHb5MmTdffddwf+n3EaWPr9aeeioiItWLBAkydP1pQpU/TUU0+publZ99xzj9OlOebs2bP6+OOPA38fP35cVVVVSklJUXZ2tpYsWaKVK1dq3LhxGjNmjJYvX67MzEzNnTvXuaIdUFhYqK1bt+r1119XcnJy4L03j8ejxMREeTweLVy4UEVFRUpJSZHb7dbixYvl9Xp1/fXXO1y9PcXFxbr55puVnZ2tpqYmbd26VW+99ZZ2797NGP1FcnJy4LMCFwwdOlSpqamB6YzTAOP0x61tePrpp012draJj483U6ZMMXv37nW6JEe9+eabRlKn24IFC4wx5y83Wr58uUlLSzMul8tMnz7dVFdXO1u0A7oaI0lmy5YtgWXOnTtn7r//fjN8+HAzZMgQc8cdd5iTJ086V7QDvve975nRo0eb+Ph4M3LkSDN9+nTzH//xH4H5jFHXPnupkTGM00DD7/kCAGBZv37PFwCAaET4AgBgGeELAIBlhC8AAJYRvgAAWEb4AgBgGeELAIBlhC8AAJYRvgAAWEb4AgBgGeELAIBl/w8k0dKrlFnRAwAAAABJRU5ErkJggg==",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.animation as animation\n",
|
|
"\n",
|
|
"fps = 60\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots(1, 1)\n",
|
|
"im = plt.imshow(numpy.linalg.norm(es[0].displacement, axis=2), cmap=\"hot\")\n",
|
|
"\n",
|
|
"def frame(i):\n",
|
|
" im.set_data(numpy.linalg.norm(es[i].displacement, axis=2))\n",
|
|
" im.autoscale()\n",
|
|
" ax.set_title(f\"$i = {i} \\qquad \\sum |F_i| = {numpy.sum(numpy.linalg.norm(es[i].force)):.3f} \\qquad \\sum |v_i| = {numpy.sum(numpy.linalg.norm(es[i].velocity)):.3f} \\qquad \\sum |d_i| = {numpy.sum(numpy.linalg.norm(es[i].displacement)):.3f}$\")\n",
|
|
" return [im]\n",
|
|
"\n",
|
|
"anim = animation.FuncAnimation(\n",
|
|
" fig,\n",
|
|
" frame, \n",
|
|
" frames = len(es),\n",
|
|
" interval = 1000 / fps,\n",
|
|
")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 42,
|
|
"id": "adc4a988-310f-4fe7-90bd-e9cd643c469e",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"anim.save('elastic.mp4', fps=fps)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "bc2f7ad3-d825-498b-84ba-7fe3b3c5eaf5",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3 (ipykernel)",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.10.11"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|