{ "cells": [ { "cell_type": "code", "execution_count": 6, "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": 7, "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": 8, "id": "b4ed52bd-02aa-411a-9af4-be9bb28227bc", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk1UlEQVR4nO3dfXCU1d3/8U9CyAMP2ZA47CY1gdQyg4gPPChEnLa35C6oQ0EZW5zYwYeBqkFFWlGq4FDFKG3VYhEqo4BTlOqMoDgVhztYrGNIIIAVsYg1P0iFDVXMLqAEzJ7fHwyrFydogps918b3a2ZnzHWd7J4cR8+Z/e7nu2nGGCMAAAAfSXc9AQAAgJNxQAEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxmuJ3A6YrGY9u7dq969eystLc31dAAAQDsYY3Tw4EEVFRUpPf3r3yNJyQPK3r17VVxc7HoaAADgNDQ2NurMM8/82jEpeUDp3bu3pON/YG5uruPZAACA9ohGoyouLo7v418nJQ8oJ8o6ubm5HFAAAEgx7fl4Bh+SBQAAvsMBBQAA+A4HFAAA4DscUAAAgO9wQAEAAL7DAQUAAPhOhw8ob7zxhsaNG6eioiKlpaVp9erVnvvGGM2ZM0eFhYXKyclReXm5du3a5Rlz4MABVVRUKDc3V3l5ebrxxht16NChb/WHAACArqPDB5TDhw/r/PPP18KFC9u8P3/+fC1YsECLFy9WbW2tevbsqTFjxujIkSPxMRUVFXr33Xe1bt06vfLKK3rjjTc0derU0/8rAABAl5JmjDGn/ctpaVq1apUmTJgg6fi7J0VFRfrVr36lX//615KkSCSiYDCoZcuWadKkSXrvvfc0aNAgbdq0ScOHD5ckrV27Vpdffrn+85//qKio6BtfNxqNKhAIKBKJdHqjttaYUV3DAe0/eER9e2drWL8+qt/9afzni0rzJYkxjGEMYxjDGM8Y16+fyL+jW3pivveuI/t3QjvJNjQ0KBwOq7y8PH4tEAhoxIgRqqmp0aRJk1RTU6O8vLz44USSysvLlZ6ertraWl155ZXW87a0tKilpSX+czQaTeS0T2nt9n2au2aH9kW+fPcnPU2KfeVIl9ejuySp+bNjjGEMYxjDGMb45vUTNaYwkK37xg3S2MGFSqaEHlDC4bAkKRgMeq4Hg8H4vXA4rL59+3onkZGh/Pz8+JiTVVVVae7cuYmc6jdau32fbv7LFp389lLspAtf/ZfIGMYwhjGMYYxfXj9RY8KRI7r5L1u06NqhST2kpESKZ9asWYpEIvFHY2Njp75ea8xo7pod1uEEAIDvmhN74dw1O9R68ommEyX0gBIKhSRJTU1NnutNTU3xe6FQSPv37/fc/+KLL3TgwIH4mJNlZWXFvxgwGV8QWNdwwFPWAQDgu8xI2hc5orqGA0l7zYQeUEpLSxUKhVRdXR2/Fo1GVVtbq7KyMklSWVmZmpubVV9fHx+zfv16xWIxjRgxIpHTOW37D3I4AQDgZMncHzv8GZRDhw7pgw8+iP/c0NCgbdu2KT8/XyUlJZo+fboeeOABDRgwQKWlpZo9e7aKioriSZ+zzz5bY8eO1ZQpU7R48WIdO3ZM06ZN06RJk9qV4EmGvr2zXU8BAADfSeb+2OEDyubNm/U///M/8Z9nzJghSZo8ebKWLVummTNn6vDhw5o6daqam5t1ySWXaO3atcrO/vKPWrFihaZNm6bRo0crPT1dEydO1IIFCxLw5yTGsH59rE82AwDwXZaednx/TJZv1QfFlc7ug1Lz7090zZKNCX9eAABS2XNTRqrsrILT/v2O7N8pkeJJNj6DAgCALZn7IweUNpzRK8v1FAAA8J1k7o8cUNqSckUvAACSIIn7IweUNnx8uOWbBwEA8B2TzP2RA0obiBkDAGBL5v7IAaUNJ2LGAADguGTHjDmgtKF+96f0QAEA4Cti5vj+mCwcUNpAzBgAABsxY8eIGQMAYCNm7BrlHQAAbMSM3SJmDACAjZixY5R4AACwUeJxjRIPAAA2SjxuUeIBAMBGiccxOskCAGCjk6xjdJIFAMCLTrI+QCdZAAC86CTrA3SSBQDARidZx4gZAwBgI2bsGuUdAABsxIzdImYMAICNmLFjxIwBALARM3aMmDEAAF7EjH2AmDEAAF7EjH2AmDEAADZixo4RMwYAwEbM2DXKOwAA2IgZu0XMGAAAGzFjxyjxAABgo8TjGiUeAABslHjcosQDAICNEo9jdJIFAMBGJ1nH6CQLAIAXnWR9gE6yAAB40UnWB+gkCwCAjU6yjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHiBkDAGAjZuwYMWMAALyIGfsAMWMAALyIGfsAMWMAAGzEjB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjlHiAQDARonHNUo8AADYKPG4RYkHAAAbJR7H6CQLAICNTrKO0UkWAAAvOsn6AJ1kAQDwopOsD9BJFgAAG51kHSNmDACAjZixa5R3AACwETN2i5gxAAA2YsaOETMGAMBGzNgxYsYAAHilfMy4tbVVs2fPVmlpqXJycnTWWWfp/vvvlzFfFq6MMZozZ44KCwuVk5Oj8vJy7dq1K9FTOW3EjAEA8Er5mPHDDz+sRYsW6U9/+pPee+89Pfzww5o/f74ef/zx+Jj58+drwYIFWrx4sWpra9WzZ0+NGTNGR474I95LzBgAAFsy98eMRD/hW2+9pfHjx+uKK66QJPXv31/PPfec6urqJB1/9+Sxxx7Tvffeq/Hjx0uSnnnmGQWDQa1evVqTJk1K9JQ6jJgxAAC2lI4ZX3zxxaqurtb7778vSXr77bf15ptv6rLLLpMkNTQ0KBwOq7y8PP47gUBAI0aMUE1NTZvP2dLSomg06nl0Kso7AADYkrg/JvwdlLvvvlvRaFQDBw5Ut27d1Nraqnnz5qmiokKSFA6HJUnBYNDze8FgMH7vZFVVVZo7d26ip3pKxIwBALCldMz4+eef14oVK/Tss89qy5YtWr58uX7/+99r+fLlp/2cs2bNUiQSiT8aGxsTOGMbJR4AAGzJ3B8T/g7KnXfeqbvvvjv+WZJzzz1Xu3fvVlVVlSZPnqxQKCRJampqUmFhYfz3mpqadMEFF7T5nFlZWcrKSuKhgRIPAAC2VO4k+9lnnyk93fu03bp1UywWkySVlpYqFAqpuro6fj8ajaq2tlZlZWWJns5pocQDAIAtmftjwt9BGTdunObNm6eSkhKdc8452rp1qx555BHdcMMNkqS0tDRNnz5dDzzwgAYMGKDS0lLNnj1bRUVFmjBhQqKnc1roJAsAgC2Z+2PCDyiPP/64Zs+erVtuuUX79+9XUVGRfvnLX2rOnDnxMTNnztThw4c1depUNTc365JLLtHatWuVne2Pg8GJTrI0awMA4Lhkd5JNM19t8ZoiotGoAoGAIpGIcnNzE/78Nf/+RNcs2Zjw5wUAIJU9N2Wkys4qOO3f78j+zXfxtIFOsgAA2JK5P3JAaQMxYwAAbCndSbZLSLmiFwAASZDKMeOugJgxAAC2lO4k2xUQMwYAwJbM/ZEDShtOxIwBAMBxyY4Zc0BpQ/3uT+mBAgDAV8TM8f0xWTigtIGYMQAANmLGjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHKPEAAGCjxOMaJR4AAGyUeNyixAMAgI0Sj2N0kgUAwEYnWcfoJAsAgBedZH2ATrIAAHjRSdYH6CQLAICNTrKOETMGAMBGzNg1yjsAANiIGbtFzBgAABsxY8eIGQMAYCNm7BgxYwAAvIgZ+wAxYwAAvIgZ+wAxYwAAbMSMHSNmDACAjZixa5R3AACwETN2i5gxAAA2YsaOUeIBAMBGicc1SjwAANgo8bhFiQcAABslHsfoJAsAgI1Oso7RSRYAAC86yfoAnWQBAPCik6wP0EkWAAAbnWQdI2YMAICNmLFrlHcAALARM3aLmDEAADZixo5R4gEAwEaJxzVKPAAA2CjxuEWJBwAAGyUex+gkCwCAjU6yjtFJFgAALzrJ+gCdZAEA8KKTrA/QSRYAABudZB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjhEzBgDARszYMWLGAAB4ETP2AWLGAAB4ETP2AWLGAADYiBk7RswYAAAbMWPXKO8AAGAjZuwWMWMAAGzEjB2jxAMAgI0Sj2uUeAAAsKV6ieejjz7Stddeq4KCAuXk5Ojcc8/V5s2b4/eNMZozZ44KCwuVk5Oj8vJy7dq1qzOmcloo8QAAYEvpEs+nn36qUaNGqXv37nr11Ve1Y8cO/eEPf1CfPl82d5k/f74WLFigxYsXq7a2Vj179tSYMWN05Ig/4r10kgUAwJbM/TEj0U/48MMPq7i4WEuXLo1fKy0tjf+zMUaPPfaY7r33Xo0fP16S9MwzzygYDGr16tWaNGlSoqfUYSc6ydKsDQCA41K+k+zLL7+s4cOH6+qrr1bfvn01ZMgQLVmyJH6/oaFB4XBY5eXl8WuBQEAjRoxQTU1Nm8/Z0tKiaDTqeXQmOskCAOCV8p1kP/zwQy1atEgDBgzQa6+9pptvvlm33Xabli9fLkkKh8OSpGAw6Pm9YDAYv3eyqqoqBQKB+KO4uDjR0/agkywAALaU7iQbi8U0dOhQPfjggxoyZIimTp2qKVOmaPHixaf9nLNmzVIkEok/GhsbEzhjGzFjAABsKR0zLiws1KBBgzzXzj77bO3Zs0eSFAqFJElNTU2eMU1NTfF7J8vKylJubq7n0ako7wAAYEvlmPGoUaO0c+dOz7X3339f/fr1k3T8A7OhUEjV1dXx+9FoVLW1tSorK0v0dE4LMWMAAGzJ3B8TnuK54447dPHFF+vBBx/Uz372M9XV1enJJ5/Uk08+KUlKS0vT9OnT9cADD2jAgAEqLS3V7NmzVVRUpAkTJiR6OqeFmDEAALaUjhlfeOGFWrVqlWbNmqXf/va3Ki0t1WOPPaaKior4mJkzZ+rw4cOaOnWqmpubdckll2jt2rXKzvbHwYCYMQAAXsmOGacZY1JuG45GowoEAopEIp3yeZSaf3+ia5ZsTPjzAgCQyp6bMlJlZxWc9u93ZP/mu3jaQMwYAABbSseMuwJixgAA2FI6ZtwlpFzRCwCAJEjlmHFXQMwYAABbSn+bcVdAiQcAABslHtco8QAAYKPE4xYlHgAAbJR4HKOTLAAAtmTujxxQ2nCikywAADgu2Z1kOaC0oX73p7S5BwDgK2Lm+P6YLBxQ2vB/O8KupwAAgO/QSdah1pjRqm0fuZ4GAAC+Q8zYobqGAzpw+JjraQAA4D/EjN3hiwIBAGgbMWOH/t/Hh11PAQAAXyJm7EhrzOi5uj2upwEAgO8QM3aoruGAwlG6yAIAcDJixg4RLwYA4NSIGTvQGjN6cSvxYgAAToWYsQN/Wr9Ln35GvBgAgFMiZpxca7fv06P/t8v1NAAA8DVixknUGjOau2aH62kAAOB7lHiSqK7hgPZFaM4GAMA3osSTPHSOBQCgfSjxJFEyu+IBAJDK6CSbRMP69VF62jePa88YAAC6KjrJJln97k8Va0dNrT1jAADoqugkm2R8BgUAgPahk2wSJTMyBQBAKiNmnEyUbgAAaB9ixsmTzMgUAACpjJhxEhEzBgCgfYgZJ9FFpfkqDGSLFDEAAKdGzDjJuqWn6b5xg/goCgAAX4OYsQNjBxfqhlH9XU8DAABfI2bswP8OCrmeAgAAvkbM2IGLSvOV3zPT9TQAAPAvYsbJ1y09TRMuKHI9DQAAfIuYsSOUeQAAODVKPI5Q5gEA4GtQ4nGDMg8AAKdGicchyjwAALSNTrIODevXR+m0lQUAwINOso7V7/5UMdrKAgDgQSdZx5LZJQ8AgFRCJ1mH+HZjAADaRszYIaLGAACcAjFjd4gaAwDQNmLGjhE1BgDARszYMaLGAAB4ETP2AaLGAAB4ETP2AaLGAADYiBk7lswYFQAAqYKYsWuUdwAAsBEzdiuZMSoAAFIFMWPHKPEAAGCjxOMaJR4AAGxdqcTz0EMPKS0tTdOnT49fO3LkiCorK1VQUKBevXpp4sSJampq6uyptBslHgAAbF2mxLNp0yb9+c9/1nnnnee5fscdd2jNmjV64YUXtGHDBu3du1dXXXVVZ06lQ/jCQAAAbF2ik+yhQ4dUUVGhJUuWqE+fLzvPRSIRPfXUU3rkkUd06aWXatiwYVq6dKneeustbdy4sbOm0yF0kgUAwKvLdJKtrKzUFVdcofLycs/1+vp6HTt2zHN94MCBKikpUU1NTZvP1dLSomg06nl0JjrJAgDglexOshmd8aQrV67Uli1btGnTJuteOBxWZmam8vLyPNeDwaDC4XCbz1dVVaW5c+d2xlTbRCdZAABsKd1JtrGxUbfffrtWrFih7OzE1KpmzZqlSCQSfzQ2NibkeU+FmDEAALaUjhnX19dr//79Gjp0qDIyMpSRkaENGzZowYIFysjIUDAY1NGjR9Xc3Oz5vaamJoVCoTafMysrS7m5uZ5Hp6K8AwCALYn7Y8JLPKNHj9Y777zjuXb99ddr4MCBuuuuu1RcXKzu3bururpaEydOlCTt3LlTe/bsUVlZWaKnc1qIGQMAYEvm/pjwA0rv3r01ePBgz7WePXuqoKAgfv3GG2/UjBkzlJ+fr9zcXN16660qKyvTyJEjEz2d00LMGAAAWzL3x075kOw3efTRR5Wenq6JEyeqpaVFY8aM0RNPPOFiKm06ETMmyQMAwHHJjhmnGWNSbhuORqMKBAKKRCKd8nmUmn9/omuW+KMnCwAAfvHclJEqO6vgtH+/I/s338XTBmLGAADYUjpm3BUQMwYAwJbSMeMuIeWKXgAAJEFX+jbjVETMGAAAW5f5NuNURYkHAAAbJR7XKPEAAGCjxOMWJR4AAGyUeByjkywAALZk7o8cUNpwopMsAAA4LtmdZDmgtKF+96e0uQcA4Cti5vj+mCwcUNpAJ1kAAGx0knWMmDEAADZixq5R3gEAwEbM2C1ixgAA2IgZO0aJBwAAGyUe1yjxAABgo8TjFiUeAABslHgco5MsAAA2Osk6RidZAAC86CTrA3SSBQDAi06yPkAnWQAAbHSSdYyYMQAANmLGrlHeAQDARszYLWLGAADYiBk7RswYAAAbMWPHiBkDAOBFzNgHiBkDAOBFzNgHiBkDAGAjZuwYMWMAAGzEjF2jvAMAgI2YsVvEjAEAsBEzdowSDwAANko8rlHiAQDARonHLUo8AADYKPE4RidZAABsdJJ1jE6yAAB40UnWB+gkCwCAF51kfYBOsgAA2Ogk6xgxYwAAbMSMXaO8AwCAjZixW8SMAQCwETN2jJgxAAA2YsaOETMGAMCLmLEPEDMGAMCLmLEPEDMGAMBGzNgxYsYAANiIGbtGeQcAABsxY7eIGQMAYCNm7BglHgAAbJR4XKPEAwCAjRKPW5R4AACwUeJxjE6yAADY6CTrGJ1kAQDwopOsD9BJFgAALzrJ+gCdZAEAsNFJ1jFixgAA2FI6ZlxVVaULL7xQvXv3Vt++fTVhwgTt3LnTM+bIkSOqrKxUQUGBevXqpYkTJ6qpqSnRUzl9lHcAALClcsx4w4YNqqys1MaNG7Vu3TodO3ZMP/nJT3T48OH4mDvuuENr1qzRCy+8oA0bNmjv3r266qqrEj2V00bMGAAAWzL3x4xEP+HatWs9Py9btkx9+/ZVfX29fvjDHyoSieipp57Ss88+q0svvVSStHTpUp199tnauHGjRo4cmegpdRgxYwAAbF0qZhyJRCRJ+fn5kqT6+nodO3ZM5eXl8TEDBw5USUmJampqOns67ULMGAAAry4VM47FYpo+fbpGjRqlwYMHS5LC4bAyMzOVl5fnGRsMBhUOh9t8npaWFkWjUc+jMxEzBgDAq0vFjCsrK7V9+3atXLnyWz1PVVWVAoFA/FFcXJygGbaNmDEAALYuETOeNm2aXnnlFb3++us688wz49dDoZCOHj2q5uZmz/impiaFQqE2n2vWrFmKRCLxR2NjY2dNWxIxYwAA2pLSMWNjjKZNm6ZVq1Zp/fr1Ki0t9dwfNmyYunfvrurq6vi1nTt3as+ePSorK2vzObOyspSbm+t5dCrKOwAA2JK4PyY8xVNZWalnn31WL730knr37h3/XEkgEFBOTo4CgYBuvPFGzZgxQ/n5+crNzdWtt96qsrIyXyR4JGLGAAC0JaVjxosWLZIk/fjHP/ZcX7p0qa677jpJ0qOPPqr09HRNnDhRLS0tGjNmjJ544olET+W0UeIBAMCWzP0x4QcUY775/Z/s7GwtXLhQCxcuTPTLJwYlHgAAbKncSbYroMQDAIAtmfsjB5Q20EkWAABbl+okm4roJAsAgFeX6iSbqugkCwCAV5fqJJuq6CQLAICtS3SSTWXEjAEAsKV0J9kugfIOAAA2YsZuETMGAMBGzNgxYsYAANiIGTtGzBgAAC9ixj5AzBgAAC9ixj5AzBgAABsxY8eIGQMAYCNm7BrlHQAAbMSM3SJmDACAjZixY5R4AACwUeJxjRIPAAA2SjxuUeIBAMBGiccxOskCAGCjk6xjdJIFAMCLTrI+QCdZAAC86CTrA3SSBQDARidZx4gZAwBgI2bsGuUdAABsxIzdImYMAICNmLFjxIwBALARM3aMmDEAAF7EjH2AmDEAAF7EjH2AmDEAADZixo4RMwYAwEbM2DXKOwAA2IgZu0XMGAAAGzFjxyjxAABgo8TjGiUeAABslHjcosQDAICNEo9jdJIFAMBGJ1nH6CQLAIAXnWR9gE6yAAB40UnWB+gkCwCAjU6yjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHiBkDAGAjZuwYMWMAALyIGfsAMWMAALyIGfsAMWMAAGzEjB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjlHiAQDARonHNUo8AADYKPG4RYkHAAAbJR7H6CQLAICNTrKO0UkWAAAvOsn6AJ1kAQDwopOsD9BJFgAAG51kHSNmDACA7TsTM164cKH69++v7OxsjRgxQnV1dS6n8yXKOwAA2L4LMeO//vWvmjFjhu677z5t2bJF559/vsaMGaP9+/e7mlIcMWMAAGzfiZjxI488oilTpuj666/XoEGDtHjxYvXo0UNPP/20qynFUeIBAMDW5Us8R48eVX19vcrLy7+cSHq6ysvLVVNTY41vaWlRNBr1PDoVJR4AAGxdvcTz8ccfq7W1VcFg0HM9GAwqHA5b46uqqhQIBOKP4uLizp0fJR4AACzfiRJPR8yaNUuRSCT+aGxs7NTXo5MsAAC2ZO6PGUl7pa8444wz1K1bNzU1NXmuNzU1KRQKWeOzsrKUlZW8utdFpfkqDGQrHDlCtQcA8J2XJikUyNZFpflJe00n76BkZmZq2LBhqq6ujl+LxWKqrq5WWVmZiyl5dEtP033jBkk6/i8FAIDvqhP74H3jBqlbEr8HxlmJZ8aMGVqyZImWL1+u9957TzfffLMOHz6s66+/3tWUPMYOLtSia4cqFPC+nXXyv5u8Ht2V16M7YxjDGMYwhjG+ev1EjQkFsrXo2qEaO7hQyeSkxCNJP//5z/Xf//5Xc+bMUTgc1gUXXKC1a9daH5x1aezgQv3voJDqGg5o/8Ej6ts7W8P69VH97k/jP594u4sxjGEMYxjDmK+Ocf36ifw7kvnOyQlpxpiU+5hFNBpVIBBQJBJRbm6u6+kAAIB26Mj+nRIpHgAA8N3CAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxxQAACA73BAAQAAvuOs1f23caL5bTQadTwTAADQXif27fY0sU/JA8rBgwclScXFxY5nAgAAOurgwYMKBAJfOyYlv4snFotp79696t27t9LSOv8LjKLRqIqLi9XY2Mh3/3Qi1jk5WOfkYJ2Tg3VOjkStszFGBw8eVFFRkdLTv/5TJin5Dkp6errOPPPMpL9ubm4u/wEkAeucHKxzcrDOycE6J0ci1vmb3jk5gQ/JAgAA3+GAAgAAfIcDSjtkZWXpvvvuU1ZWluupdGmsc3KwzsnBOicH65wcLtY5JT8kCwAAujbeQQEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcED5BgsXLlT//v2VnZ2tESNGqK6uzvWUUlpVVZUuvPBC9e7dW3379tWECRO0c+dOz5gjR46osrJSBQUF6tWrlyZOnKimpiZHM+4aHnroIaWlpWn69Onxa6xzYnz00Ue69tprVVBQoJycHJ177rnavHlz/L4xRnPmzFFhYaFycnJUXl6uXbt2OZxx6mltbdXs2bNVWlqqnJwcnXXWWbr//vs93+fCOnfcG2+8oXHjxqmoqEhpaWlavXq153571vTAgQOqqKhQbm6u8vLydOONN+rQoUOJmaDBKa1cudJkZmaap59+2rz77rtmypQpJi8vzzQ1NbmeWsoaM2aMWbp0qdm+fbvZtm2bufzyy01JSYk5dOhQfMxNN91kiouLTXV1tdm8ebMZOXKkufjiix3OOrXV1dWZ/v37m/POO8/cfvvt8eus87d34MAB069fP3PdddeZ2tpa8+GHH5rXXnvNfPDBB/ExDz30kAkEAmb16tXm7bffNj/96U9NaWmp+fzzzx3OPLXMmzfPFBQUmFdeecU0NDSYF154wfTq1cv88Y9/jI9hnTvub3/7m7nnnnvMiy++aCSZVatWee63Z03Hjh1rzj//fLNx40bzj3/8w/zgBz8w11xzTULmxwHla1x00UWmsrIy/nNra6spKioyVVVVDmfVtezfv99IMhs2bDDGGNPc3Gy6d+9uXnjhhfiY9957z0gyNTU1rqaZsg4ePGgGDBhg1q1bZ370ox/FDyisc2Lcdddd5pJLLjnl/VgsZkKhkPnd734Xv9bc3GyysrLMc889l4wpdglXXHGFueGGGzzXrrrqKlNRUWGMYZ0T4eQDSnvWdMeOHUaS2bRpU3zMq6++atLS0sxHH330redEiecUjh49qvr6epWXl8evpaenq7y8XDU1NQ5n1rVEIhFJUn5+viSpvr5ex44d86z7wIEDVVJSwrqfhsrKSl1xxRWe9ZRY50R5+eWXNXz4cF199dXq27evhgwZoiVLlsTvNzQ0KBwOe9Y5EAhoxIgRrHMHXHzxxaqurtb7778vSXr77bf15ptv6rLLLpPEOneG9qxpTU2N8vLyNHz48PiY8vJypaenq7a29lvPISW/LDAZPv74Y7W2tioYDHquB4NB/etf/3I0q64lFotp+vTpGjVqlAYPHixJCofDyszMVF5enmdsMBhUOBx2MMvUtXLlSm3ZskWbNm2y7rHOifHhhx9q0aJFmjFjhn7zm99o06ZNuu2225SZmanJkyfH17Kt/4+wzu139913KxqNauDAgerWrZtaW1s1b948VVRUSBLr3Anas6bhcFh9+/b13M/IyFB+fn5C1p0DCpyprKzU9u3b9eabb7qeSpfT2Nio22+/XevWrVN2drbr6XRZsVhMw4cP14MPPihJGjJkiLZv367Fixdr8uTJjmfXdTz//PNasWKFnn32WZ1zzjnatm2bpk+frqKiIta5C6PEcwpnnHGGunXrZqUampqaFAqFHM2q65g2bZpeeeUVvf766zrzzDPj10OhkI4eParm5mbPeNa9Y+rr67V//34NHTpUGRkZysjI0IYNG7RgwQJlZGQoGAyyzglQWFioQYMGea6dffbZ2rNnjyTF15L/j3w7d955p+6++25NmjRJ5557rn7xi1/ojjvuUFVVlSTWuTO0Z01DoZD279/vuf/FF1/owIEDCVl3DiinkJmZqWHDhqm6ujp+LRaLqbq6WmVlZQ5nltqMMZo2bZpWrVql9evXq7S01HN/2LBh6t69u2fdd+7cqT179rDuHTB69Gi988472rZtW/wxfPhwVVRUxP+Zdf72Ro0aZcXk33//ffXr10+SVFpaqlAo5FnnaDSq2tpa1rkDPvvsM6Wne7erbt26KRaLSWKdO0N71rSsrEzNzc2qr6+Pj1m/fr1isZhGjBjx7SfxrT9m24WtXLnSZGVlmWXLlpkdO3aYqVOnmry8PBMOh11PLWXdfPPNJhAImL///e9m37598cdnn30WH3PTTTeZkpISs379erN582ZTVlZmysrKHM66a/hqiscY1jkR6urqTEZGhpk3b57ZtWuXWbFihenRo4f5y1/+Eh/z0EMPmby8PPPSSy+Zf/7zn2b8+PHEXzto8uTJ5nvf+148Zvziiy+aM844w8ycOTM+hnXuuIMHD5qtW7earVu3GknmkUceMVu3bjW7d+82xrRvTceOHWuGDBliamtrzZtvvmkGDBhAzDhZHn/8cVNSUmIyMzPNRRddZDZu3Oh6SilNUpuPpUuXxsd8/vnn5pZbbjF9+vQxPXr0MFdeeaXZt2+fu0l3EScfUFjnxFizZo0ZPHiwycrKMgMHDjRPPvmk534sFjOzZ882wWDQZGVlmdGjR5udO3c6mm1qikaj5vbbbzclJSUmOzvbfP/73zf33HOPaWlpiY9hnTvu9ddfb/P/x5MnTzbGtG9NP/nkE3PNNdeYXr16mdzcXHP99debgwcPJmR+acZ8pRUfAACAD/AZFAAA4DscUAAAgO9wQAEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxxQAACA73BAAQAAvsMBBQAA+M7/B3qz1s2r8+4UAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "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(100, 100)\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": null, "id": "0f1f9dfd-8266-4874-9ba5-83a8a784d0e9", "metadata": {}, "outputs": [], "source": [ "import multiprocessing\n", "\n", "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", "\n", " def job_force(row):\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", " with multiprocessing.Pool(4) as pool:\n", " pool.map(new_force, range(elm.n_rows))\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": 9, "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": 10, "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.1))" ] }, { "cell_type": "code", "execution_count": 12, "id": "421c0fd2-2dba-4fd5-89ff-0968227b08b9", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "7d7f0752567340c08c97df22788fc638", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=0, description='i', max=1000), Output()), _dom_classes=('widget-interact…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ " None>" ] }, "execution_count": 12, "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": null, "id": "96ab8cf1-8739-4faa-b6b5-70aa52a15a2a", "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.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }