NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4fa62fbd-0920-4990-892e-29ef20062ee1",
   "metadata": {},
   "source": [
    "# Elastodynamics with Newmark time-stepping\n",
    "\n",
    "$u$ ... displacement, $v = \\dot u $ ... velocity, $a = \\dot v$ .. acceleration\n",
    "\n",
    "Newmark scheme:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{u^{n+1}-u^n}{\\tau} = \\frac{v^n+v^{n+1}}{2} \\\\\n",
    "\\frac{v^{n+1}-v^n}{\\tau} = \\frac{a^n+a^{n+1}}{2} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "with new acceleration, with elasticity operator $K$:\n",
    "\n",
    "$$\n",
    "a^{n+1} = f - K(u^{n+1})\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9e13aaa-6c05-4759-b4f0-b4f25a685439",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e35649b7-f3a2-4641-bf33-ab3d3aeb6e2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2707f0b34351481eacd732441f9a7607",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from netgen.occ import *\n",
    "shape = Rectangle(1,0.1).Face()\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "shape.edges.Min(Y).name=\"bot\"\n",
    "shape.vertices.Min(X+Y).maxh=0.01\n",
    "shape.vertices.Min(X-Y).maxh=0.01\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a225848b-237b-47a5-ac94-a3d4f05f24e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = 210, 0.2\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "def C(u):\n",
    "    F = Id(2) + Grad(u)\n",
    "    return F.trans * F\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)\n",
    "\n",
    "tau = 0.02\n",
    "tend = 20\n",
    "force = CF( (0, -1) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "17ef5205-2520-4739-8814-69b744f5ee17",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = VectorH1(mesh, order=3, dirichlet=\"left\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfv = GridFunction(fes)\n",
    "gfa = GridFunction(fes)\n",
    "gfuold = GridFunction(fes)\n",
    "gfvold = GridFunction(fes)\n",
    "gfaold = GridFunction(fes)\n",
    "\n",
    "bfa = BilinearForm(fes)\n",
    "bfa += Variation(NeoHooke(C(u))*dx)\n",
    "\n",
    "vel_new = 2/tau * (u-gfuold) - gfvold\n",
    "acc_new = 2/tau * (vel_new-gfvold) - gfaold\n",
    "\n",
    "bfa += acc_new*v*dx\n",
    "bfa += -force*v*dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "76b97cb7-3f75-4efb-8c5a-6bb308cf32b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6b8e07f3a0b14dc5b88cd8f8e4ebb5e1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ab746565e1046968bf1a1f77c16a6ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gfu_history = GridFunction(fes, multidim=0)\n",
    "sceneu = Draw (gfu, deformation=True)\n",
    "scenev = Draw (gfv)\n",
    "gfu.vec[:] = 0\n",
    "t = 0\n",
    "step = 0\n",
    "while t < tend:\n",
    "    t += tau\n",
    "    step += 1\n",
    "    solvers.Newton(a=bfa, u=gfu, printing=False, inverse=\"sparsecholesky\")\n",
    "\n",
    "    gfv.vec[:] = 2/tau * (gfu.vec-gfuold.vec) - gfvold.vec\n",
    "    gfa.vec[:] = 2/tau * (gfv.vec-gfvold.vec) - gfaold.vec\n",
    "\n",
    "    sceneu.Redraw()\n",
    "    scenev.Redraw()\n",
    "    \n",
    "    gfuold.vec[:] = gfu.vec\n",
    "    gfvold.vec[:] = gfv.vec\n",
    "    gfaold.vec[:] = gfa.vec\n",
    "    if step%20 == 0: gfu_history.AddMultiDimComponent(gfu.vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "51060e8b-2688-4a18-bd5f-5da141c3b35c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b2defdb94cb44b30950cf4827349b649",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Draw (gfu_history, mesh, interpolate_multidim=True, animate=True, min=0, max=1, autoscale=False, deformation=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e275e2a0-e8ee-4f47-948f-92235c38f967",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}