Paul Brinkmeier
251ae05fd7
Currently buggy, there seems to be a mix-up in the ordering of x and y coordinates somewhere. This probably stems from the use of row+column ordering.
397 lines
30 KiB
Plaintext
397 lines
30 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 179,
|
|
"id": "e465f900-0a94-46d4-b3c5-23cc0129557b",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy\n",
|
|
"\n",
|
|
"from typing import Iterator, Tuple"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7e2231fc-95a1-427a-a557-3b01e4fec19a",
|
|
"metadata": {},
|
|
"source": [
|
|
"(O'Brien, 2014)\n",
|
|
"\n",
|
|
"$E_i \\in \\mathbb{R}$ is the elastic energy for node $i$:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"E_i\n",
|
|
"= \\sum^N_{j = 1} \\frac{K_{ij}}{2} (\\mathbf{u}_{ij} \\cdot \\mathbf{n}_{ij})^2\n",
|
|
"+ \\sum^{N^2}_{jk} \\frac{c_{ijk}}{2} \\left( \\cos \\theta_{jik} - \\frac{1}{\\sqrt{2}} \\right)^2\n",
|
|
"$$\n",
|
|
"\n",
|
|
"where\n",
|
|
"\n",
|
|
"* $K_{ij} \\in \\mathbb{R}$ is the elastic spring constant\n",
|
|
"* $N = 8$ is the number of neighbors (for the 2D lattice)\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",
|
|
"* $\\mathbf{n}_{ij} = \\frac{\\mathbf{x}_{ij}}{|\\mathbf{x}_{ij}|} \\in \\mathbb{R}^2$\n",
|
|
"* $c_{jik} \\in \\mathbb{R}$ is the bond-bending constant\n",
|
|
"* $\\theta_{jik} \\in \\mathbb{R}$ is the angle between $j$ and $k$ with $i$ as the apex\n",
|
|
"\n",
|
|
"The first sum is over all the neighbors $j$ of $i$.\n",
|
|
"The second sum is over *all* pairs of neighbors $jk$ of $i$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 309,
|
|
"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 __init__(self, n_rows: int, n_cols: int):\n",
|
|
" self.n_rows = n_rows\n",
|
|
" self.n_cols = n_cols\n",
|
|
" self.displacement = numpy.zeros((n_rows, n_cols, 2))\n",
|
|
" self.velocity = numpy.zeros((n_rows, n_cols, 2))\n",
|
|
" self.force = numpy.zeros((n_rows, n_cols, 2))\n",
|
|
"\n",
|
|
" self.elastic_spring_constant = 1.0\n",
|
|
" self.bond_bending_constant = 1.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",
|
|
" # 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}\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}$ 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 Lamé constants are given by\n",
|
|
"\n",
|
|
"$\\lambda_\\text{2D} = K - \\frac{c}{\\Delta x^2}$\n",
|
|
"\n",
|
|
"and\n",
|
|
"\n",
|
|
"$\\mu_\\text{2D} = K + \\frac{c}{\\Delta x^2}$\n",
|
|
"\n",
|
|
"where $\\Delta x$ is the lattice grid spacing."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 310,
|
|
"id": "b4ed52bd-02aa-411a-9af4-be9bb28227bc",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([-2., 0.])"
|
|
]
|
|
},
|
|
"execution_count": 310,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"def neighbors(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",
|
|
" if row_offset == 0 and col_offset == 0:\n",
|
|
" continue\n",
|
|
" else:\n",
|
|
" k = 1.0 if row_offset == 0 or col_offset == 0 else 2.0\n",
|
|
" yield (k, row + row_offset, col + col_offset)\n",
|
|
"\n",
|
|
"def calculate_force(elm: ELM, 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 _, nb_row, nb_col in neighbors(row, col):\n",
|
|
" u = elm.get_displacement(row, col) - elm.get_displacement(nb_row, nb_col)\n",
|
|
" x = elm.get_lattice_position(row, col) - elm.get_lattice_position(nb_row, nb_col)\n",
|
|
"\n",
|
|
" # print((k, nb_row, nb_col))\n",
|
|
" \n",
|
|
" total -= K * numpy.dot(u, x) + (c * u) / numpy.linalg.norm(x) ** 2\n",
|
|
"\n",
|
|
" return total\n",
|
|
"\n",
|
|
"e = ELM(20, 20)\n",
|
|
"for row in range(1, e.n_rows - 1):\n",
|
|
" e.set_displacement(row, 2, numpy.array([0.5, 0.0]))\n",
|
|
"calculate_force(e, 2, 2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 311,
|
|
"id": "94d7da0c-a01a-4040-95e0-8f8b913d273a",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"(2.0, 1, 1)\n",
|
|
"(1.0, 1, 2)\n",
|
|
"(2.0, 1, 3)\n",
|
|
"(1.0, 2, 1)\n",
|
|
"(1.0, 2, 3)\n",
|
|
"(2.0, 3, 1)\n",
|
|
"(1.0, 3, 2)\n",
|
|
"(2.0, 3, 3)\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1t0lEQVR4nO3de1hVdd7//9cGZaMGWxnllISYRqOmlAdCM7VQMCOZ+7qn9Fep3dpMjnaNY4eJ7gybE2pOZZODTWNR06jllHqPFWYaOibq7YFbzfKbhmfAyYm9gRQV1u8Px5VbDrKRw0d8Pq5rXbLXeq+1Px8XH9aLtdZeOCzLsgQAAGAwv+ZuAAAAwKUQWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxmvV3A1oCJWVlTp27JiCgoLkcDiauzkAAKAOLMtSSUmJIiMj5edX+zmUFhFYjh07pqioqOZuBgAAqIfDhw+rc+fOtda0iMASFBQk6VyHg4ODm7k1AACgLjwej6KiouzjeG1aRGA5fxkoODiYwAIAwBWmLrdzcNMtAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGC8FvHguMZSUWlpS/6/dLzklEKDAjUgJkT+fvytIqCpMAaB5mfKOPQpsGRmZiozM1MHDhyQJPXs2VPPPvusRo4cWeM6S5cu1YwZM3TgwAF1795ds2fP1l133WUvtyxL6enpeu2111RcXKxBgwYpMzNT3bt3r1+PGkj27gI99/c9KnCfsudFuAKVntJDyb0imrFlwNWBMQg0P5PGoU+XhDp37qxZs2Zp27Zt2rp1q+644w6NHj1an3/+ebX1Gzdu1NixYzVx4kTt2LFDqampSk1N1e7du+2aOXPm6OWXX9aCBQu0efNmtWvXTklJSTp16lS122wK2bsLNPnt7V47SJIK3ac0+e3tyt5d0EwtA64OjEGg+Zk2Dh2WZVmXs4GQkBA9//zzmjhxYpVl9913n8rKyrRy5Up73q233qq4uDgtWLBAlmUpMjJSjz32mB5//HFJktvtVlhYmLKysjRmzJg6tcHj8cjlcsntdl/23xKqqLR02+y19g466y6SZVlq3T5ckuSQFO4K1IZf3sGpaaARXDwGK0+VquJUqfzbBMvP2ZYxCDSBC8ehZVk66y6Sw89PrYJDJTXcsdCX43e9b7qtqKjQkiVLVFZWpoSEhGprcnNzlZiY6DUvKSlJubm5kqT8/HwVFhZ61bhcLsXHx9s11SkvL5fH4/GaGsqW/H/ZPyhPfr1NRxdM1LFXJ+m7/3euPZakAvcpbcn/V4O9J4DvXTgGLatShYvTdOzVSTr+t5nn5okxCDS2C8dh2e61OvbqJB3NnKjyo19Iap5x6HNg2bVrl6655ho5nU498sgjWrZsmXr06FFtbWFhocLCwrzmhYWFqbCw0F5+fl5NNdXJyMiQy+Wyp6ioKF+7UaPjJd+f+irZvvKCrz+osQ5Aw7lwbJUf2aMzx/Ptr08f/7raOgANy+tYuOP8sdBSyY4Pa6xrbD4HltjYWOXl5Wnz5s2aPHmyxo8frz179jRG22qUlpYmt9ttT4cPH26wbYcGBdpfB/f/0QVfp9ZYB6DhXDi2nJ17KiDi3A34gdG9FRDatdo6AA3L61jYL/XfXzkU1G90jXWNzeePNQcEBKhbt26SpL59++p///d/NW/ePL366qtVasPDw1VUVOQ1r6ioSOHh4fby8/MiIiK8auLi4mpsg9PplNPp9LXpdTIgJkQRrkAVuk8pMLq3Oj/6V8my5N+uvaTvr9sNiAlplPcHrnYXjkE5HAobkyHrzEk5AtpKYgwCTeHCcdiuxxAFRveWHH7yb+uS1Dzj8LIfHFdZWany8vJqlyUkJGjNmjVe81avXm3f8xITE6Pw8HCvGo/Ho82bN9d4X0xj8/dzKD3l3CUuhyT/ti6vsCJJ6Sk9uNkPaCQXj0G/gED5t+sgv9ZOxiDQRKocC9t18AorUtOPQ58CS1pamtavX68DBw5o165dSktLU05Oju6//35J0rhx45SWlmbX//znP1d2drZ+//vf68svv9TMmTO1detWTZ06VZLkcDg0bdo0/eY3v9H//M//aNeuXRo3bpwiIyOVmpracL30UXKvCGU+cIvCXd6nusJdgcp84BaeAQE0MsYg0PxMG4c+XRI6fvy4xo0bp4KCArlcLvXu3VurVq3S8OHDJUmHDh2Sn9/3GWjgwIFatGiRnnnmGT399NPq3r27li9frl69etk1Tz75pMrKyvSTn/xExcXFuu2225Sdna3AwOa9Pp3cK0LDe4Qb8XQ/4GrEGASan0nj8LKfw2KChnwOCwAAaBpN8hwWAACApkJgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwnk+BJSMjQ/3791dQUJBCQ0OVmpqqvXv31rrO0KFD5XA4qkyjRo2yayZMmFBleXJycv16BAAAWpxWvhSvW7dOU6ZMUf/+/XX27Fk9/fTTGjFihPbs2aN27dpVu87777+v06dP269PnDihPn366Mc//rFXXXJyst544w37tdPp9KVpAACgBfMpsGRnZ3u9zsrKUmhoqLZt26bbb7+92nVCQkK8Xi9ZskRt27atElicTqfCw8N9aQ4AALhKXNY9LG63W1LVUFKbhQsXasyYMVXOyOTk5Cg0NFSxsbGaPHmyTpw4UeM2ysvL5fF4vCYAANByOSzLsuqzYmVlpe655x4VFxdrw4YNdVpny5Ytio+P1+bNmzVgwAB7/vmzLjExMdq/f7+efvppXXPNNcrNzZW/v3+V7cycOVPPPfdclflut1vBwcH16Q4AAGhiHo9HLperTsfvegeWyZMn66OPPtKGDRvUuXPnOq3z05/+VLm5udq5c2etdV9//bWuv/56ffLJJ7rzzjurLC8vL1d5ebn92uPxKCoqisACAMAVxJfAUq9LQlOnTtXKlSv16aef1jmslJWVacmSJZo4ceIla7t27aqOHTtq37591S53Op0KDg72mgAAQMvl0023lmXp0Ucf1bJly5STk6OYmJg6r7t06VKVl5frgQceuGTtkSNHdOLECUVERPjSPAAA0EL5dIZlypQpevvtt7Vo0SIFBQWpsLBQhYWFOnnypF0zbtw4paWlVVl34cKFSk1N1Q9+8AOv+aWlpXriiSe0adMmHThwQGvWrNHo0aPVrVs3JSUl1bNbAACgJfHpDEtmZqakcw+Du9Abb7yhCRMmSJIOHTokPz/vHLR3715t2LBBH3/8cZVt+vv7a+fOnXrzzTdVXFysyMhIjRgxQr/+9a95FgsAAJB0GTfdmsSXm3YAAIAZGv2mWwAAgKZEYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMJ5PgSUjI0P9+/dXUFCQQkNDlZqaqr1799a6TlZWlhwOh9cUGBjoVWNZlp599llFRESoTZs2SkxM1FdffeV7bwAAQIvkU2BZt26dpkyZok2bNmn16tU6c+aMRowYobKyslrXCw4OVkFBgT0dPHjQa/mcOXP08ssva8GCBdq8ebPatWunpKQknTp1yvceAQCAFqeVL8XZ2dler7OyshQaGqpt27bp9ttvr3E9h8Oh8PDwapdZlqWXXnpJzzzzjEaPHi1JeuuttxQWFqbly5drzJgxvjQRAAC0QJd1D4vb7ZYkhYSE1FpXWlqq6OhoRUVFafTo0fr888/tZfn5+SosLFRiYqI9z+VyKT4+Xrm5udVur7y8XB6Px2sCAAAtV70DS2VlpaZNm6ZBgwapV69eNdbFxsbq9ddf14oVK/T222+rsrJSAwcO1JEjRyRJhYWFkqSwsDCv9cLCwuxlF8vIyJDL5bKnqKio+nYDAABcAeodWKZMmaLdu3dryZIltdYlJCRo3LhxiouL05AhQ/T++++rU6dOevXVV+v71kpLS5Pb7banw4cP13tbAADAfD7dw3Le1KlTtXLlSq1fv16dO3f2ad3WrVvr5ptv1r59+yTJvrelqKhIERERdl1RUZHi4uKq3YbT6ZTT6axP0wEAwBXIpzMslmVp6tSpWrZsmdauXauYmBif37CiokK7du2yw0lMTIzCw8O1Zs0au8bj8Wjz5s1KSEjwefsAAKDl8ekMy5QpU7Ro0SKtWLFCQUFB9j0mLpdLbdq0kSSNGzdO1157rTIyMiRJv/rVr3TrrbeqW7duKi4u1vPPP6+DBw9q0qRJks59gmjatGn6zW9+o+7duysmJkYzZsxQZGSkUlNTG7CrAADgSuVTYMnMzJQkDR061Gv+G2+8oQkTJkiSDh06JD+/70/cfPvtt3r44YdVWFioDh06qG/fvtq4caN69Ohh1zz55JMqKyvTT37yExUXF+u2225TdnZ2lQfMAQCAq5PDsiyruRtxuTwej1wul9xut4KDg5u7OQAAoA58OX7zt4QAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMbzKbBkZGSof//+CgoKUmhoqFJTU7V3795a13nttdc0ePBgdejQQR06dFBiYqK2bNniVTNhwgQ5HA6vKTk52ffeAACAFsmnwLJu3TpNmTJFmzZt0urVq3XmzBmNGDFCZWVlNa6Tk5OjsWPH6tNPP1Vubq6ioqI0YsQIHT161KsuOTlZBQUF9rR48eL69QgAALQ4DsuyrPqu/M9//lOhoaFat26dbr/99jqtU1FRoQ4dOuiVV17RuHHjJJ07w1JcXKzly5fXqx0ej0cul0tut1vBwcH12gYAAGhavhy/L+seFrfbLUkKCQmp8zrfffedzpw5U2WdnJwchYaGKjY2VpMnT9aJEydq3EZ5ebk8Ho/XBAAAWq56n2GprKzUPffco+LiYm3YsKHO6/3sZz/TqlWr9PnnnyswMFCStGTJErVt21YxMTHav3+/nn76aV1zzTXKzc2Vv79/lW3MnDlTzz33XJX5nGEBAODK4csZlnoHlsmTJ+ujjz7Shg0b1Llz5zqtM2vWLM2ZM0c5OTnq3bt3jXVff/21rr/+en3yySe68847qywvLy9XeXm5/drj8SgqKorAAgDAFaTRLwlNnTpVK1eu1KefflrnsDJ37lzNmjVLH3/8ca1hRZK6du2qjh07at++fdUudzqdCg4O9poAAEDL1cqXYsuy9Oijj2rZsmXKyclRTExMndabM2eOfvvb32rVqlXq16/fJeuPHDmiEydOKCIiwpfmAQCAFsqnMyxTpkzR22+/rUWLFikoKEiFhYUqLCzUyZMn7Zpx48YpLS3Nfj179mzNmDFDr7/+urp06WKvU1paKkkqLS3VE088oU2bNunAgQNas2aNRo8erW7duikpKamBugkAAK5kPgWWzMxMud1uDR06VBEREfb0zjvv2DWHDh1SQUGB1zqnT5/Wf/7nf3qtM3fuXEmSv7+/du7cqXvuuUc33HCDJk6cqL59++of//iHnE5nA3UTAABcyS7rOSym4DksAABceZrsOSwAAABNgcACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCyX4ezZs83dBOCqdvbsWVmW1dzNAK5qTXUsJLDUoqLSUu7+E1qRd1S5+0+ootL7B+NLL72k5ORkvfLKKzp48GAztRJowSorpPx/SLv+du7fygqvxcXFxerbt6+mT5+utWvX6syZM83UUKDlutSx8JlnntHo0aP15z//WQUFBY3WDoflw68nGRkZev/99/Xll1+qTZs2GjhwoGbPnq3Y2Nha11u6dKlmzJihAwcOqHv37po9e7buuusue7llWUpPT9drr72m4uJiDRo0SJmZmerevXud2uXxeORyueR2uxUcHFzX7tQqe3eBnvv7HhW4T9nzIlyBSk/poeReEfb7dunSRd9++60k6aabblJKSopSUlLUv39/+fv7N0hbgKvSnv+Rsn8peY59Py84UkqeLfW4x56VlpamWbNmSZJcLpeSk5OVkpKi5ORk/eAHP2jqVgMtSl2OhUePHlXXrl11+vRpSVK/fv3sY2FcXJwcDkeN2/fl+O1TYElOTtaYMWPUv39/nT17Vk8//bR2796tPXv2qF27dtWus3HjRt1+++3KyMjQ3XffrUWLFmn27Nnavn27evXqJUmaPXu2MjIy9OabbyomJkYzZszQrl27tGfPHgUGBl6yXQ0dWLJ3F2jy29t18X/M+f/yzAdusXfUb3/7Wz3zzDNVttGpUyeNGjVKKSkpGj58uIKCgi67XcBVY8//SO+Ok2oahfe+ZYeWb775Rl26dFFZWZlXpZ+fnwYNGmT/4IyNja31BycAb74cCx999FG98sorVbZx7bXX6u6771ZKSoruuOMOtWnTxmt5owWWi/3zn/9UaGio1q1bp9tvv73amvvuu09lZWVauXKlPe/WW29VXFycFixYIMuyFBkZqccee0yPP/64JMntdissLExZWVkaM2bMJdvRkIGlotLSbbPXeqXJ0p2rdeEPTleb1kob+UP5+Tnk8Xj0i1/8otZtBgQEaOjQofYPzujo6MtqI9CiVVZIL/XyOrOy9ViFdhadvxzkkALbS8N/Jfmdu6q9cOFCbdy4sdbNduvWzf7BOXjwYLVu3bqROgBc+S4+FlpnT6tsT45XzYXHwqNHj+rZZ5+tdZtt2rTR8OHDdffdd+vuu+9WRERE0wWWffv2qXv37tq1a5d9tuRi1113naZPn65p06bZ89LT07V8+XL93//9n77++mtdf/312rFjh+Li4uyaIUOGKC4uTvPmzauyzfLycpWXl9uvPR6PoqKiGiSw5O4/obGvbfKad/D50VWunV+O85eO7r77bg0YMIBLR8CF8v8hvXm316y0T05p1menG+wtzl86uvvuuzVy5EguHQEXufhYWHHSoyMv/38N+h79+vXT8OHDlZGRUafjd71vuq2srNS0adM0aNCgGsOKJBUWFiosLMxrXlhYmAoLC+3l5+fVVHOxjIwMuVwue4qKiqpvN6o4XnLq0kWXadeuXfrd736ngQMHatCgQdq6dWujvydwxSgtavS3cLvdeuedd/Tggw8qOjpac+bMsa+/A2iaY+HWrVuVkZFR5/pW9X2jKVOmaPfu3dqwYUN9N1FvaWlpmj59uv36/BmWhhAadOl7Zi6Hw+FQQkKCfWq6Z8+eXFcHLnRN2KVrLlNYWJg9BhMTE2u8Bw+4WjX2sdDf31+DBw/W8OHD9d///d91WqdegWXq1KlauXKl1q9fr86dO9daGx4erqIi79+YioqKFB4ebi8/Py8iIsKr5sJLRBdyOp1yOp31afolDYgJUYQrUIXuU/ZdKx1TnpD+feXMIal929b61ehe8vNz6Ntvv9UjjzxS6zavueYaJSUlKSUlRXfddZc6derUKG0HWoTogec+DeQp0Pl7x8be1Fo3R1xw6bRNiHTXXPselldeeUX/+Mc/at1sXFycfR9Z37595efHUx2Amlx8LPRr3UYd7/mlvfziY+HBgwf15JNP1rrNDh06aOTIkUpJSVFSUpI6dOggj8fTOIHFsiw9+uijWrZsmXJychQTE3PJdRISErRmzRqve1hWr16thIQESVJMTIzCw8O1Zs0aO6B4PB5t3rxZkydP9qV5DcLfz6H0lB6a/PZ2OXTux2W7G2+T9P2d0S9fcGd0enp6tdvp0qWLfZ/KkCFDGi1gAS2On/+5jy6/O0769yjsHeav3mH++v5TQn+0PyVUWFioCRMmVNmM0+nUnXfead/g15CXjoGW7uJjoVq1VrsfDpZU/bHwpz/9abXbiY2Ntc9mDho0SK1a1fvCjm833f7sZz/TokWLtGLFCq9nr7hcLvujSuPGjdO1115rX5fauHGjhgwZolmzZmnUqFFasmSJfve731X5WPOsWbO8Pta8c+fOZvtYs1S3z55/++236tKlizwejxwOh2699Vb7Nzgu9QCXqdrnsFwrJc/yeg7LY489phdeeEESl3qAhlaXY+HBgwfVrVs3nT171r7Uc/5YeKnnqTXap4RqOgC/8cYb9m84Q4cOVZcuXZSVlWUvX7p0qZ555hn7wXFz5syp9sFxf/rTn1RcXKzbbrtNf/zjH3XDDTfUqV2NEVikcx/r2pL/Lx0vOaXQoEANiAmRv9/3/wdz587Vpk2blJKSopEjRyo0NLTB3huAzn067+DGczfiXhN27nKR3/eXhoqKivQf//Ef9pmUfv36cakHaGCXOhampaXp0KFDXpd66qrJPtZsisYKLJdSWVnJD0egGTEGgeZ3OePQl+M3I/0y8IMSaF6MQaD5NdU4ZLQDAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABjP58Cyfv16paSkKDIyUg6HQ8uXL6+1fsKECXI4HFWmnj172jUzZ86ssvzGG2/0uTMAAKBl8jmwlJWVqU+fPpo/f36d6ufNm6eCggJ7Onz4sEJCQvTjH//Yq65nz55edRs2bPC1aQAAoIVq5esKI0eO1MiRI+tc73K55HK57NfLly/Xt99+q4ceesi7Ia1aKTw83NfmAACAq0CT38OycOFCJSYmKjo62mv+V199pcjISHXt2lX333+/Dh06VOM2ysvL5fF4vCYAANByNWlgOXbsmD766CNNmjTJa358fLyysrKUnZ2tzMxM5efna/DgwSopKal2OxkZGfaZG5fLpaioqKZoPgAAaCYOy7Kseq/scGjZsmVKTU2tU31GRoZ+//vf69ixYwoICKixrri4WNHR0XrhhRc0ceLEKsvLy8tVXl5uv/Z4PIqKipLb7VZwcLDP/QAAAE3P4/HI5XLV6fjt8z0s9WVZll5//XU9+OCDtYYVSWrfvr1uuOEG7du3r9rlTqdTTqezMZoJAAAM1GSXhNatW6d9+/ZVe8bkYqWlpdq/f78iIiKaoGUAAMB0PgeW0tJS5eXlKS8vT5KUn5+vvLw8+ybZtLQ0jRs3rsp6CxcuVHx8vHr16lVl2eOPP65169bpwIED2rhxo370ox/J399fY8eO9bV5AACgBfL5ktDWrVs1bNgw+/X06dMlSePHj1dWVpYKCgqqfMLH7Xbrvffe07x586rd5pEjRzR27FidOHFCnTp10m233aZNmzapU6dOvjYPAAC0QJd1060pfLlpBwAAmMGX4zd/SwgAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDyfA8v69euVkpKiyMhIORwOLV++vNb6nJwcORyOKlNhYaFX3fz589WlSxcFBgYqPj5eW7Zs8bVpAACghfI5sJSVlalPnz6aP3++T+vt3btXBQUF9hQaGmove+eddzR9+nSlp6dr+/bt6tOnj5KSknT8+HFfmwcAAFqgVr6uMHLkSI0cOdLnNwoNDVX79u2rXfbCCy/o4Ycf1kMPPSRJWrBggT744AO9/vrreuqpp3x+LwAA0LI02T0scXFxioiI0PDhw/XZZ5/Z80+fPq1t27YpMTHx+0b5+SkxMVG5ubnVbqu8vFwej8drAgAALVejB5aIiAgtWLBA7733nt577z1FRUVp6NCh2r59uyTpm2++UUVFhcLCwrzWCwsLq3Kfy3kZGRlyuVz2FBUV1djdAAAAzcjnS0K+io2NVWxsrP164MCB2r9/v1588UX95S9/qdc209LSNH36dPu1x+MhtAAA0II1emCpzoABA7RhwwZJUseOHeXv76+ioiKvmqKiIoWHh1e7vtPplNPpbPR2AgAAMzTLc1jy8vIUEREhSQoICFDfvn21Zs0ae3llZaXWrFmjhISE5mgeAAAwjM9nWEpLS7Vv3z77dX5+vvLy8hQSEqLrrrtOaWlpOnr0qN566y1J0ksvvaSYmBj17NlTp06d0p///GetXbtWH3/8sb2N6dOna/z48erXr58GDBigl156SWVlZfanhgAAwNXN58CydetWDRs2zH59/l6S8ePHKysrSwUFBTp06JC9/PTp03rsscd09OhRtW3bVr1799Ynn3zitY377rtP//znP/Xss8+qsLBQcXFxys7OrnIjLgAAuDo5LMuymrsRl8vj8cjlcsntdis4OLi5mwMAAOrAl+M3f0sIAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8nwPL+vXrlZKSosjISDkcDi1fvrzW+vfff1/Dhw9Xp06dFBwcrISEBK1atcqrZubMmXI4HF7TjTfe6GvTAABAC+VzYCkrK1OfPn00f/78OtWvX79ew4cP14cffqht27Zp2LBhSklJ0Y4dO7zqevbsqYKCAnvasGGDr00DAAAtVCtfVxg5cqRGjhxZ5/qXXnrJ6/Xvfvc7rVixQn//+9918803f9+QVq0UHh7ua3MAAMBVoMnvYamsrFRJSYlCQkK85n/11VeKjIxU165ddf/99+vQoUM1bqO8vFwej8drAgAALVeTB5a5c+eqtLRU9957rz0vPj5eWVlZys7OVmZmpvLz8zV48GCVlJRUu42MjAy5XC57ioqKaqrmAwCAZuCwLMuq98oOh5YtW6bU1NQ61S9atEgPP/ywVqxYocTExBrriouLFR0drRdeeEETJ06ssry8vFzl5eX2a4/Ho6ioKLndbgUHB/vcDwAA0PQ8Ho9cLledjt8+38NSX0uWLNGkSZO0dOnSWsOKJLVv31433HCD9u3bV+1yp9Mpp9PZGM0EAAAGapJLQosXL9ZDDz2kxYsXa9SoUZesLy0t1f79+xUREdEErQMAAKbz+QxLaWmp15mP/Px85eXlKSQkRNddd53S0tJ09OhRvfXWW5LOXQYaP3685s2bp/j4eBUWFkqS2rRpI5fLJUl6/PHHlZKSoujoaB07dkzp6eny9/fX2LFjG6KPAADgCufzGZatW7fq5ptvtj+SPH36dN1888169tlnJUkFBQVen/D505/+pLNnz2rKlCmKiIiwp5///Od2zZEjRzR27FjFxsbq3nvv1Q9+8ANt2rRJnTp1utz+AQCAFuCybro1hS837QAAADP4cvzmbwkBAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABivVXM3wGQVlZa25P9Lx0tOKTQoUANiQuTv52juZgFXDcYg0PxMGYc+B5b169fr+eef17Zt21RQUKBly5YpNTW11nVycnI0ffp0ff7554qKitIzzzyjCRMmeNXMnz9fzz//vAoLC9WnTx/94Q9/0IABA3xtXoPJ3l2g5/6+RwXuU/a8CFeg0lN6KLlXRLO1C7haMAaB5mfSOPT5klBZWZn69Omj+fPn16k+Pz9fo0aN0rBhw5SXl6dp06Zp0qRJWrVqlV3zzjvvaPr06UpPT9f27dvVp08fJSUl6fjx4742r0Fk7y7Q5Le3e+0gSSp0n9Lkt7cre3dBs7QLuFowBoHmZ9o4dFiWZdV7ZYfjkmdYfvnLX+qDDz7Q7t277XljxoxRcXGxsrOzJUnx8fHq37+/XnnlFUlSZWWloqKi9Oijj+qpp566ZDs8Ho9cLpfcbreCg4Pr2x1J50593TZ7rb2DzrqLZFmWWrcPlyQ5JIW7ArXhl3dwahpoBBePwcpTpao4VSr/NsHyc7ZlDAJN4MJxaFmWzrqL5PDzU6vgUEkNdyz05fjd6Dfd5ubmKjEx0WteUlKScnNzJUmnT5/Wtm3bvGr8/PyUmJho11ysvLxcHo/Ha2ooW/L/Zf+gPPn1Nh1dMFHHXp2k7/7fubZYkgrcp7Ql/18N9p4AvnfhGLSsShUuTtOxVyfp+N9mnpsnxiDQ2C4ch2W71+rYq5N0NHOiyo9+Ial5xmGjB5bCwkKFhYV5zQsLC5PH49HJkyf1zTffqKKiotqawsLCareZkZEhl8tlT1FRUQ3W3uMl35/6Ktm+8oKvP6ixDkDDuXBslR/ZozPH8+2vTx//uto6AA3L61i44/yx0FLJjg9rrGtsV+THmtPS0uR2u+3p8OHDDbbt0KBA++vg/j+64OvUGusANJwLx5azc08FRHSXJAVG91ZAaNdq6wA0LK9jYb/Uf3/lUFC/0TXWNbZG/1hzeHi4ioqKvOYVFRUpODhYbdq0kb+/v/z9/autCQ8Pr3abTqdTTqezUdo7ICZEEa5AFbpPKTC6tzo/+lfJsuTfrr2k76/bDYgJaZT3B652F45BORwKG5Mh68xJOQLaSmIMAk3hwnHYrscQBUb3lhx+8m/rktQ847DRz7AkJCRozZo1XvNWr16thIQESVJAQID69u3rVVNZWak1a9bYNU3J38+h9JQeks7tEP+2Lq+wIknpKT242Q9oJBePQb+AQPm36yC/1k7GINBEqhwL23XwCitS049DnwNLaWmp8vLylJeXJ+ncx5bz8vJ06NAhSecu14wbN86uf+SRR/T111/rySef1Jdffqk//vGPevfdd/WLX/zCrpk+fbpee+01vfnmm/riiy80efJklZWV6aGHHrrM7tVPcq8IZT5wi8Jd3qe6wl2BynzgFp4BATQyxiDQ/Ewbhz5/rDknJ0fDhg2rMn/8+PHKysrShAkTdODAAeXk5Hit84tf/EJ79uxR586dNWPGjCoPjnvllVfsB8fFxcXp5ZdfVnx8fJ3a1JAfa76QKU/3A65WjEGg+TXmOPTl+H1Zz2ExRWMFFgAA0HiMeg4LAADA5SKwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGa/S/1twUzj+s1+PxNHNLAABAXZ0/btflofstIrCUlJRIkqKiopq5JQAAwFclJSVyuVy11rSIvyVUWVmpY8eOKSgoSA5Hw/5hNI/Ho6ioKB0+fLhF/p2ilt4/qeX3kf5d+Vp6H1t6/6SW38fG6p9lWSopKVFkZKT8/Gq/S6VFnGHx8/NT586dG/U9goODW+Q34XktvX9Sy+8j/bvytfQ+tvT+SS2/j43Rv0udWTmPm24BAIDxCCwAAMB4BJZLcDqdSk9Pl9PpbO6mNIqW3j+p5feR/l35WnofW3r/pJbfRxP61yJuugUAAC0bZ1gAAIDxCCwAAMB4BBYAAGA8AgsAADDeVRVY1q9fr5SUFEVGRsrhcGj58uWXXCcnJ0e33HKLnE6nunXrpqysrCo18+fPV5cuXRQYGKj4+Hht2bKl4RtfR7728f3339fw4cPVqVMnBQcHKyEhQatWrfKqmTlzphwOh9d04403NmIvauZr/3Jycqq03eFwqLCw0KvOlH3oa/8mTJhQbf969uxp15i0/zIyMtS/f38FBQUpNDRUqamp2rt37yXXW7p0qW688UYFBgbqpptu0ocffui13LIsPfvss4qIiFCbNm2UmJior776qrG6Uav69PG1117T4MGD1aFDB3Xo0EGJiYlVvger29fJycmN2ZVq1ad/WVlZVdoeGBjoVWPKPqxP/4YOHVrtOBw1apRdY8r+k6TMzEz17t3bfghcQkKCPvroo1rXMWEMXlWBpaysTH369NH8+fPrVJ+fn69Ro0Zp2LBhysvL07Rp0zRp0iSvA/o777yj6dOnKz09Xdu3b1efPn2UlJSk48ePN1Y3auVrH9evX6/hw4frww8/1LZt2zRs2DClpKRox44dXnU9e/ZUQUGBPW3YsKExmn9JvvbvvL1793q1PzQ01F5m0j70tX/z5s3z6tfhw4cVEhKiH//4x151puy/devWacqUKdq0aZNWr16tM2fOaMSIESorK6txnY0bN2rs2LGaOHGiduzYodTUVKWmpmr37t12zZw5c/Tyyy9rwYIF2rx5s9q1a6ekpCSdOnWqKbrlpT59zMnJ0dixY/Xpp58qNzdXUVFRGjFihI4ePepVl5yc7LUfFy9e3NjdqaI+/ZPOPSH1wrYfPHjQa7kp+7A+/Xv//fe9+rZ79275+/tXGYcm7D9J6ty5s2bNmqVt27Zp69atuuOOOzR69Gh9/vnn1dYbMwatq5Qka9myZbXWPPnkk1bPnj295t13331WUlKS/XrAgAHWlClT7NcVFRVWZGSklZGR0aDtrY+69LE6PXr0sJ577jn7dXp6utWnT5+Ga1gDqUv/Pv30U0uS9e2339ZYY+o+rM/+W7ZsmeVwOKwDBw7Y80zdf5ZlWcePH7ckWevWraux5t5777VGjRrlNS8+Pt766U9/almWZVVWVlrh4eHW888/by8vLi62nE6ntXjx4sZpuA/q0seLnT171goKCrLefPNNe9748eOt0aNHN0ILL09d+vfGG29YLperxuUm78P67L8XX3zRCgoKskpLS+15pu6/8zp06GD9+c9/rnaZKWPwqjrD4qvc3FwlJiZ6zUtKSlJubq4k6fTp09q2bZtXjZ+fnxITE+2aK01lZaVKSkoUEhLiNf+rr75SZGSkunbtqvvvv1+HDh1qphbWT1xcnCIiIjR8+HB99tln9vyWtg8XLlyoxMRERUdHe803df+53W5JqvL9dqFLjcP8/HwVFhZ61bhcLsXHxxuxD+vSx4t99913OnPmTJV1cnJyFBoaqtjYWE2ePFknTpxo0LbWR137V1paqujoaEVFRVX5bd7kfVif/bdw4UKNGTNG7dq185pv4v6rqKjQkiVLVFZWpoSEhGprTBmDBJZaFBYWKiwszGteWFiYPB6PTp48qW+++UYVFRXV1lx8j8SVYu7cuSotLdW9995rz4uPj1dWVpays7OVmZmp/Px8DR48WCUlJc3Y0rqJiIjQggUL9N577+m9995TVFSUhg4dqu3bt0tSi9qHx44d00cffaRJkyZ5zTd1/1VWVmratGkaNGiQevXqVWNdTePw/P45/6+J+7CufbzYL3/5S0VGRnodAJKTk/XWW29pzZo1mj17ttatW6eRI0eqoqKiMZpeJ3XtX2xsrF5//XWtWLFCb7/9tiorKzVw4EAdOXJEkrn7sD77b8uWLdq9e3eVcWja/tu1a5euueYaOZ1OPfLII1q2bJl69OhRba0pY7BF/LVmNIxFixbpueee04oVK7zu8Rg5cqT9de/evRUfH6/o6Gi9++67mjhxYnM0tc5iY2MVGxtrvx44cKD279+vF198UX/5y1+asWUN780331T79u2VmprqNd/U/TdlyhTt3r272e6naQr16eOsWbO0ZMkS5eTkeN2YOmbMGPvrm266Sb1799b111+vnJwc3XnnnQ3a7rqqa/8SEhK8fnsfOHCgfvjDH+rVV1/Vr3/968ZuZr3VZ/8tXLhQN910kwYMGOA137T9Fxsbq7y8PLndbv3tb3/T+PHjtW7duhpDiwk4w1KL8PBwFRUVec0rKipScHCw2rRpo44dO8rf37/amvDw8KZs6mVbsmSJJk2apHfffbfKqb+LtW/fXjfccIP27dvXRK1rWAMGDLDb3lL2oWVZev311/Xggw8qICCg1loT9t/UqVO1cuVKffrpp+rcuXOttTWNw/P75/y/pu1DX/p43ty5czVr1ix9/PHH6t27d621Xbt2VceOHZttP9anf+e1bt1aN998s912E/dhffpXVlamJUuW1OkXgebefwEBAerWrZv69u2rjIwM9enTR/Pmzau21pQxSGCpRUJCgtasWeM1b/Xq1fZvCgEBAerbt69XTWVlpdasWVPjtUATLV68WA899JAWL17s9TG8mpSWlmr//v2KiIhogtY1vLy8PLvtLWUfrlu3Tvv27avTD8rm3H+WZWnq1KlatmyZ1q5dq5iYmEuuc6lxGBMTo/DwcK8aj8ejzZs3N8s+rE8fpXOfsvj1r3+t7Oxs9evX75L1R44c0YkTJ5p8P9a3fxeqqKjQrl277LabtA8vp39Lly5VeXm5HnjggUvWNtf+q0llZaXKy8urXWbMGGyw23evACUlJdaOHTusHTt2WJKsF154wdqxY4d18OBBy7Is66mnnrIefPBBu/7rr7+22rZtaz3xxBPWF198Yc2fP9/y9/e3srOz7ZolS5ZYTqfTysrKsvbs2WP95Cc/sdq3b28VFhY2ef8sy/c+/vWvf7VatWplzZ8/3yooKLCn4uJiu+axxx6zcnJyrPz8fOuzzz6zEhMTrY4dO1rHjx83vn8vvviitXz5cuurr76ydu3aZf385z+3/Pz8rE8++cSuMWkf+tq/8x544AErPj6+2m2atP8mT55suVwuKycnx+v77bvvvrNrHnzwQeupp56yX3/22WdWq1atrLlz51pffPGFlZ6ebrVu3dratWuXXTNr1iyrffv21ooVK6ydO3dao0ePtmJiYqyTJ082af8sq359nDVrlhUQEGD97W9/81qnpKTEsqxz3xePP/64lZuba+Xn51uffPKJdcstt1jdu3e3Tp06ZXz/nnvuOWvVqlXW/v37rW3btlljxoyxAgMDrc8//9yuMWUf1qd/5912223WfffdV2W+SfvPss79HFm3bp2Vn59v7dy503rqqacsh8Nhffzxx5ZlmTsGr6rAcv4jrhdP48ePtyzr3MfOhgwZUmWduLg4KyAgwOratav1xhtvVNnuH/7wB+u6666zAgICrAEDBlibNm1q/M7UwNc+DhkypNZ6yzr3Ue6IiAgrICDAuvbaa6377rvP2rdvX9N27N987d/s2bOt66+/3goMDLRCQkKsoUOHWmvXrq2yXVP2YX2+R4uLi602bdpYf/rTn6rdpkn7r7q+SfIaV0OGDPH6/rMsy3r33XetG264wQoICLB69uxpffDBB17LKysrrRkzZlhhYWGW0+m07rzzTmvv3r1N0KOq6tPH6OjoatdJT0+3LMuyvvvuO2vEiBFWp06drNatW1vR0dHWww8/3Cyhuj79mzZtmj2+wsLCrLvuusvavn2713ZN2Yf1/R798ssvLUn2Qf9CJu0/y7Ks//qv/7Kio6OtgIAAq1OnTtadd97p1W5Tx6DDsiyrgU7WAAAANAruYQEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeP8/gDXt29hNAtUAAAAASUVORK5CYII=",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"nbs = numpy.array([(ro, co, e.get_displacement(ro, co)[0], e.get_displacement(ro, co)[1]) for (k, ro, co) in neighbors(2, 2)])\n",
|
|
"\n",
|
|
"plt.scatter(nbs[:, 0], nbs[:, 1])\n",
|
|
"plt.scatter([2], [2])\n",
|
|
"plt.quiver(nbs[:, 0], nbs[:, 1], nbs[:, 2], nbs[:, 3], scale=10)\n",
|
|
"plt.quiver([2], [2], e.get_displacement(2, 2)[0], e.get_displacement(2, 2)[1], scale=10)\n",
|
|
"\n",
|
|
"for x in neighbors(2, 2):\n",
|
|
" print(x)"
|
|
]
|
|
},
|
|
{
|
|
"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": 312,
|
|
"id": "ff7adc30-191e-4b7c-b146-2a533463d6ea",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def velocity_verlet_step(elm: ELM, delta_t: float) -> ELM:\n",
|
|
" new_elm = ELM(elm.n_rows, elm.n_cols)\n",
|
|
"\n",
|
|
" new_elm.displacement = elm.displacement + elm.velocity * delta_t + 0.5 * elm.force * delta_t * delta_t\n",
|
|
" for row in range(1, elm.n_rows - 1):\n",
|
|
" for col in range(1, elm.n_cols - 1):\n",
|
|
" new_elm.set_force(row, col, calculate_force(new_elm, row, col))\n",
|
|
" new_elm.velocity = elm.velocity + 0.5 * (elm.force + new_elm.force) * delta_t\n",
|
|
"\n",
|
|
" return new_elm"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 313,
|
|
"id": "bea964bf-6e94-4b72-a458-50c2964bbba1",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"es = [e]\n",
|
|
"for _ in range(1000):\n",
|
|
" es.append(velocity_verlet_step(es[-1], 0.002))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 314,
|
|
"id": "ceb8c029-6697-41eb-9581-2a5901d7716a",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"%matplotlib inline"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 315,
|
|
"id": "674b00a2-6d40-41e6-9b0d-fd12b4847daa",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from ipywidgets import interact"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 316,
|
|
"id": "421c0fd2-2dba-4fd5-89ff-0968227b08b9",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "029ce9d43a78401786d7cec5dd94cb5c",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"interactive(children=(IntSlider(value=500, description='i', max=1000), Output()), _dom_classes=('widget-intera…"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<function __main__.do_plot(i: int) -> None>"
|
|
]
|
|
},
|
|
"execution_count": 316,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"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",
|
|
"def do_plot(i: int) -> None:\n",
|
|
" plt.scatter(\n",
|
|
" x + es[i].displacement[:, :, 0],\n",
|
|
" y + es[i].displacement[:, :, 1]\n",
|
|
" )\n",
|
|
" plt.quiver(\n",
|
|
" x + es[i].displacement[:, :, 0],\n",
|
|
" y + es[i].displacement[:, :, 1],\n",
|
|
" es[i].force[:, :, 0],\n",
|
|
" es[i].force[:, :, 1],\n",
|
|
" scale = 100\n",
|
|
" )\n",
|
|
" plt.show()\n",
|
|
"\n",
|
|
"interact(do_plot, i = (0, len(es) - 1, 1))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "64b10510-de4f-4ff6-a38f-3cb1db71a444",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "0c448b53-6820-4c7b-844b-379c21822792",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "a8fdbabe-1aee-4579-ad24-8a3284fa0372",
|
|
"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
|
|
}
|