{ "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": [ "
" ] }, "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": [ " 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 }