From 2d88aa1976d202a056e8221833dd6986f5451ad7 Mon Sep 17 00:00:00 2001 From: Paul Brinkmeier Date: Sun, 18 Jun 2023 19:10:18 +0200 Subject: [PATCH] Add notebook for debugging ELM force terms --- experiments/ELM_Debug.ipynb | 242 ++++++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) create mode 100644 experiments/ELM_Debug.ipynb diff --git a/experiments/ELM_Debug.ipynb b/experiments/ELM_Debug.ipynb new file mode 100644 index 0000000..b255e5b --- /dev/null +++ b/experiments/ELM_Debug.ipynb @@ -0,0 +1,242 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2342d5ae-1839-4132-ace8-da8449223a1c", + "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": 33, + "id": "05d4e22c-3d7d-4ac2-bf22-59b78aa40b84", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from ipywidgets import interact\n", + "from typing import Tuple" + ] + }, + { + "cell_type": "markdown", + "id": "50ba2e84-bb40-44c8-bb86-8804fd26f674", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{F}_{ij} = K_{ij} (\\mathbf{u}_{ij} \\cdot \\mathbf{x}_{ij}) + \\frac{c \\mathbf{u}_{ij}}{|\\mathbf{x}_{ij}|^2}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "id": "82bb4a02-1ecf-4c84-8b9d-b753dee75337", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "j=(0, 0) K_ij=2.0000000000000004 u_ij=array([0.1, 0. ]) x_ij=array([1., 1.])\n" + ] + }, + { + "data": { + "text/plain": [ + "-0.20000000000000007" + ] + }, + "execution_count": 61, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def calculate_force_elastic(\n", + " displacement: np.ndarray,\n", + " i: Tuple[int, int],\n", + " j: Tuple[int, int]\n", + ") -> float:\n", + " u_ij = displacement[i] - displacement[j]\n", + " x_ij = np.array(i, dtype=float) - np.array(j, dtype=float)\n", + " K_ij = np.linalg.norm(x_ij) ** 2\n", + "\n", + " print(f\"{j=} {K_ij=} {u_ij=} {x_ij=}\")\n", + "\n", + " return -K_ij * np.dot(u_ij, x_ij)\n", + "\n", + "def calculate_force_bond(\n", + " displacement: np.ndarray,\n", + " i: Tuple[int, int],\n", + " j: Tuple[int, int]\n", + ") -> float:\n", + " return 0.0\n", + "\n", + "def calculate_force(\n", + " displacement: np.ndarray,\n", + " i: Tuple[int, int],\n", + " j: Tuple[int, int]\n", + ") -> float:\n", + " return calculate_force_elastic(displacement, i, j) + calculate_force_bond(displacement, i, j)\n", + "\n", + "displacement = np.array([\n", + " [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]],\n", + " [[0.0, 0.0], [0.1, 0.0], [0.0, 0.0]],\n", + " [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]\n", + "])\n", + "calculate_force_elastic(displacement, (1, 1), (0, 0))" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "f7a019c2-8d20-4e84-9e0c-ef788781a07a", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "89618a87c5514c6698d62e3344092fdf", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='dp_x', max=1.0, min=-1.0, step=0.05), FloatSlider(va…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " None>" + ] + }, + "execution_count": 71, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def do_plot(dp_x: float, dp_y: float) -> None:\n", + " i = (1, 1)\n", + " displacement[i] = np.array([dp_x, dp_y])\n", + " \n", + " width, height, _ = displacement.shape\n", + "\n", + " xs = []\n", + " ys = []\n", + " elastic_xs = []\n", + " elastic_ys = []\n", + " \n", + " for x in range(width):\n", + " for y in range(height):\n", + " j = (x, y)\n", + " xs.append(x + displacement[j][0])\n", + " ys.append(y + displacement[j][1])\n", + "\n", + " if j != i:\n", + " F_ij = calculate_force_elastic(displacement, i, j)\n", + " F_ij_vec = (\n", + " (i + displacement[i])\n", + " - (j + displacement[j])\n", + " ) * F_ij\n", + " print(f\"{F_ij_vec=}\")\n", + " elastic_xs.append(F_ij_vec[0])\n", + " elastic_ys.append(F_ij_vec[1])\n", + "\n", + " pi = i + displacement[i]\n", + " plt.scatter(xs, ys)\n", + " plt.quiver(\n", + " [pi[0]] * len(elastic_xs),\n", + " [pi[1]] * len(elastic_ys),\n", + " elastic_xs,\n", + " elastic_ys,\n", + " scale = 5,\n", + " pivot = \"tip\",\n", + " width = 0.005\n", + " )\n", + " plt.quiver(\n", + " pi[0], pi[1],\n", + " sum(elastic_xs),\n", + " sum(elastic_ys),\n", + " pivot = \"tip\",\n", + " scale = 5,\n", + " color = \"red\"\n", + " )\n", + " \n", + "\n", + "interact(\n", + " do_plot,\n", + " dp_x = (-1.0, 1.0, 0.05),\n", + " dp_y = (-1.0, 1.0, 0.05)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4c55b68-903a-48de-8226-9cfbbbca4af6", + "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 +}