NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe42a41b-64a0-4a2f-95ce-200d3d3dfc43",
   "metadata": {},
   "source": [
    "# Shells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb492468-82b7-4f40-a602-de55a62d6425",
   "metadata": {},
   "source": [
    "Shells are the “primadonna of structures” (Ramm, 1986). On the one hand, if properly designed, they can exhibit an optimal ratio of stiffness to weight. On the other hand, especially for thin shells, they can be extremely sensitive to imperfections and develope unstable behavior.\n",
    "\n",
    "In this chapter we present efficient discretization methods to simulate nonlinear shells.\n",
    "\n",
    "* [Introduction of shell energies and (non-)linear shell models](linear_nonlinear_shell_equ.ipynb)\n",
    "* [Linear shell models](linear_KL_RM_shell_HHJ_TDNNS.ipynb)\n",
    "* [Nonlinear shell models](nonlinear_Koiter_Naghdi_shell.ipynb)\n",
    "* [Shell buckling](torus_buckling.ipynb)\n",
    "\n",
    "__Example: Sydney opera hall__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b03230d0-d1e8-4804-a969-50f43202d7ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58af7c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first fanel\n",
    "v1l = Pnt(0, -10, 0.1)\n",
    "v1t = Pnt(-6, 0, 10.5)\n",
    "w1 = Wire(ArcOfCircle(v1l, Vec(-0.3, 0.4, 1), v1t))\n",
    "\n",
    "v12l = Pnt(0.25, -10, 0.1)\n",
    "v12t = Pnt(-0.5, 0, 11)\n",
    "w12 = Wire(ArcOfCircle(v12l, Vec(-0.25, 0.3, 1), v12t))\n",
    "\n",
    "v2l = Pnt(0.5, -10, 0.1)\n",
    "v2t = Pnt(5 - 0.15, 0, 10)\n",
    "w2 = Wire(ArcOfCircle(v2l, Vec(0, 0.2, 1), v2t))\n",
    "\n",
    "# vent\n",
    "v1lv = Pnt(1, -10.5, 1.5)\n",
    "v1tv = Pnt(5, 0, 10)\n",
    "w2v = Wire(ArcOfCircle(v1lv, Vec(0, 0.4, 1), v1tv))\n",
    "\n",
    "v2lv = Pnt(4.5, -10.5, 0.5)\n",
    "v2tv = v1tv\n",
    "w3v = Wire(ArcOfCircle(v2lv, Vec(0, 0.2, 1), v2tv))\n",
    "\n",
    "v3l = Pnt(5, -10, 1)\n",
    "v3t = Pnt(5 + 0.15, 0, 10)\n",
    "w3 = Wire(ArcOfCircle(v3l, Vec(0, 0, 1), v3t))\n",
    "\n",
    "v4l = Pnt(5.5, -10 - 1, 1)\n",
    "v4t = Pnt(3.5, 0, 15)\n",
    "w4 = Wire(ArcOfCircle(v4l, Vec(0, 0, 1), v4t))\n",
    "\n",
    "# second panel\n",
    "v212l = Pnt(5.5 + 0.25, -10 - 1, 1)\n",
    "v212t = Pnt(3.5 + 5.5, 0, 11 + 2)\n",
    "w212 = Wire(ArcOfCircle(v212l, Vec(0, 0, 1), v212t))\n",
    "\n",
    "v22l = Pnt(5.5 + 0.5, -10 - 1, 1)\n",
    "v22t = Pnt(3.5 + 11 - 0.15, 0, 10 + 1)\n",
    "w22 = Wire(ArcOfCircle(v22l, Vec(0.2, 0, 1), v22t))\n",
    "\n",
    "# vent\n",
    "v21lv = Pnt(5.5 + 0.5 + 0.5, -10 - 1.5, 1.5)\n",
    "v21tv = Pnt(3.5 + 11, 0, 10 + 1)\n",
    "w22v = Wire(ArcOfCircle(v21lv, Vec(0.3, 0.1, 1), v21tv))\n",
    "\n",
    "v22lv = Pnt(3.5 + 11 - 1.5, -10 - 1.5, 1)\n",
    "v22tv = v21tv\n",
    "w23v = Wire(ArcOfCircle(v22lv, Vec(-0.2, 0, 1), v22tv))\n",
    "\n",
    "\n",
    "v23l = Pnt(3.5 + 11, -10 - 1, 0)\n",
    "v23t = Pnt(3.5 + 11 + 0.15, 0, 10 + 1)\n",
    "w23 = Wire(ArcOfCircle(v23l, Vec(0, 0, 1), v23t))\n",
    "\n",
    "v24l = Pnt(3.5 + 11, -10 - 1.5, 0)\n",
    "v24t = Pnt(11, 0, 14 + 8)\n",
    "w24 = Wire(ArcOfCircle(v24l, Vec(-0.1, 0, 1), v24t))\n",
    "\n",
    "# third panel\n",
    "v312l = Pnt(3.5 + 11 + 0.25, -10 - 1.5, 0)\n",
    "v312t = Pnt(11 + 9, 0, 11 + 6.5)\n",
    "w312 = Wire(ArcOfCircle(v312l, Vec(0.2, 0, 1), v312t))\n",
    "\n",
    "v32l = Pnt(3.5 + 11 + 0.5, -10 - 1.5, 0)\n",
    "v32t = Pnt(11 + 15 - 0.25, 0, 9)\n",
    "w32 = Wire(ArcOfCircle(v32l, Vec(1, 0, 1), v32t))\n",
    "\n",
    "# fouth panel\n",
    "\n",
    "# vent\n",
    "v32lv = Pnt(3.5 + 11 + 0.5 + 4, -10 - 1.5 - 0.5, 3)\n",
    "v32tv = Pnt(11 + 15 - 0.25, 0, 9)\n",
    "w32v = Wire(ArcOfCircle(v32lv, Vec(0.6, 0.3, 1), v32tv))\n",
    "\n",
    "v33lv = Pnt(3.5 + 11 + 0.5 + 9, -10 - 1.5 - 0.5, 2)\n",
    "v33tv = Pnt(11 + 15, 0, 9)\n",
    "w33v = Wire(ArcOfCircle(v33lv, Vec(0, 0.3, 1), v33tv))\n",
    "\n",
    "v34lv = Pnt(3.5 + 11 + 17 - 2, -10 - 1.5 - 0.5, 3)\n",
    "v34tv = Pnt(11 + 15, 0, 9)\n",
    "w34v = Wire(ArcOfCircle(v34lv, Vec(-0.6, 0.3, 1), v34tv))\n",
    "\n",
    "\n",
    "v41l = Pnt(3.5 + 11 + 17, -10 - 1.5, 0)\n",
    "v41t = Pnt(11 + 15 + 0.25, 0, 9)\n",
    "w41 = Wire(ArcOfCircle(v41l, Vec(-0.5, 0, 1), v41t))\n",
    "\n",
    "v42l = Pnt(3.5 + 11 + 17 + 0.25, -10 - 1.5, 0)\n",
    "v42t = Pnt(3.5 + 11 + 17 + 0.25 - 0.5, 0, 11.75)\n",
    "w42 = Wire(ArcOfCircle(v42l, Vec(0, 0, 1), v42t))\n",
    "\n",
    "v43l = Pnt(3.5 + 11 + 17 + 0.5, -10 - 1.5, 0)\n",
    "v43t = Pnt(3.5 + 11 + 17 + 0.5 + 5, 0, 12)\n",
    "w43 = Wire(ArcOfCircle(v43l, Vec(0.25, 0, 1), v43t))\n",
    "\n",
    "shape = Glue(\n",
    "    [w1, w12, w2, w2v, w3v, w3, w4]\n",
    "    + [w212, w22, w22v, w23v, w23, w24]\n",
    "    + [w312, w32, w32v, w33v, w34v]\n",
    "    + [w41, w42, w43]\n",
    ")\n",
    "Draw(shape);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c2ec76f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = ThruSections([w2, w12, w1], solid=False)\n",
    "f2 = ThruSections([w2v, w2], solid=False)\n",
    "v1 = ThruSections([w3v, w2v], solid=False)\n",
    "v2 = ThruSections([w3, w3v], solid=False)\n",
    "f3 = ThruSections([w4, w3], solid=False)\n",
    "\n",
    "f4 = ThruSections([w22, w212, w4], solid=False)\n",
    "f5 = ThruSections([w22v, w22], solid=False)\n",
    "v3 = ThruSections([w23v, w22v], solid=False)\n",
    "v4 = ThruSections([w23, w23v], solid=False)\n",
    "f6 = ThruSections([w24, w23], solid=False)\n",
    "\n",
    "f7 = ThruSections([w32, w312, w24], solid=False)\n",
    "v5 = ThruSections([w32v, w32], solid=False)\n",
    "v6 = ThruSections([w33v, w32v], solid=False)\n",
    "v7 = ThruSections([w34v, w33v], solid=False)\n",
    "\n",
    "f8 = ThruSections([w41, w34v], solid=False)\n",
    "f9 = ThruSections([w43, w42, w41], solid=False)\n",
    "shell = Glue([f1, f2, f3, f4, f5, f6, f7, f8, f9, v1, v2, v3, v4, v5, v6, v7])\n",
    "\n",
    "shell = Glue([shell, shell.Mirror(Axes((0, 0, 0), Y))])\n",
    "\n",
    "shell.edges[Z <= 0.2].name = \"fix\"\n",
    "\n",
    "shell.vertices.hpref = 1\n",
    "shell.edges.hpref = 1\n",
    "shell.faces.col = (222/255, 224/255, 204/255,1)\n",
    "\n",
    "Draw(shell)\n",
    "print (\"(model created by Sebastian Hirnschall)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54939ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(shell).GenerateMesh(maxh=2))\n",
    "mesh.RefineHP(1)\n",
    "mesh.Curve(3)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adef764c",
   "metadata": {},
   "outputs": [],
   "source": [
    "thickness = 0.1\n",
    "order = 4\n",
    "\n",
    "E = 2.85e4\n",
    "nu = 0.3\n",
    "mu = E / (2 * (1 + nu))\n",
    "lam = E * nu / (1 - nu**2)\n",
    "\n",
    "\n",
    "def MaterialNorm(mat, E, nu):\n",
    "    return E / (1 - nu**2) * ((1 - nu) * InnerProduct(mat, mat) + nu * Trace(mat) ** 2)\n",
    "\n",
    "\n",
    "def MaterialNormInv(mat, E, nu):\n",
    "    return (1 + nu) / E * (InnerProduct(mat, mat) - nu / (nu + 1) * Trace(mat) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "120f420f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesU = VectorH1(mesh, order=order, dirichlet_bbnd=\"fix\")\n",
    "fesS = HDivDivSurface(mesh, order=order - 1, discontinuous=True)\n",
    "fesH = NormalFacetSurface(mesh, order=order - 1)  # , dirichlet_bbnd=\"fix\")\n",
    "\n",
    "fes = fesU * fesS * fesH\n",
    "(u, sigma, uh) = fes.TrialFunction()\n",
    "# Need to take traces as we are on the surface\n",
    "sigma, uh = sigma.Trace(), uh.Trace()\n",
    "\n",
    "fesRegge = HCurlCurl(mesh, order=order - 1)\n",
    "\n",
    "# normal, tangential, and co-normal vectors\n",
    "nv = specialcf.normal(3)\n",
    "tv = specialcf.tangential(3)\n",
    "cnv = Cross(nv, tv)\n",
    "\n",
    "# Projection to tangent space and surface derivatives\n",
    "Ptau = Id(3) - OuterProduct(nv, nv)\n",
    "grad_u = Grad(u).Trace()\n",
    "Ftau = grad_u + Ptau\n",
    "Etautau = 0.5 * (Ftau.trans * Ftau - Ptau)\n",
    "\n",
    "# sum_{i=1}^3 hesse(u_i) nv_i\n",
    "Hn = (u.Operator(\"hesseboundary\").trans * nv).Reshape((3, 3))\n",
    "\n",
    "bfa = BilinearForm(fes, symmetric=True, condense=True)\n",
    "# membrane part\n",
    "bfa += Variation(\n",
    "    thickness / 2 * MaterialNorm(Interpolate(Etautau, fesRegge), E, nu) * ds\n",
    ")\n",
    "# bending part\n",
    "bfa += Variation(-12 / thickness**3 / 2 * MaterialNormInv(sigma, E, nu) * ds)\n",
    "bfa += Variation(InnerProduct(sigma, Hn) * ds)\n",
    "bfa += Variation(\n",
    "    sigma[cnv, cnv] * (uh - grad_u[nv, :]) * cnv * ds(element_boundary=True)\n",
    ")\n",
    "\n",
    "force = CF((0, 0, -0.01))\n",
    "bfa += Variation(-force * u * ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ac3913",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "\n",
    "with TaskManager():\n",
    "    solvers.Newton(a=bfa, u=gfu, maxerr=1e-8, printing=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f338bf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(gfu.components[0], scale=10, deformation=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0152fa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(Norm(gfu.components[1]), mesh, order=3, max=0.2);"
   ]
  }
 ],
 "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
}