{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.1.4 Shell model\n",
    "## Simple Naghdi shell model\n",
    "Geometric model and meshing. Clamped on left boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import *\n",
    "from ngsolve import *\n",
    "from ngsolve.internal import visoptions\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "order = 3\n",
    "\n",
    "geo = CSGeometry()\n",
    "cyl   = Cylinder(Pnt(0,0,0),Pnt(1,0,0),0.4).bc(\"cyl\")\n",
    "left  = Plane(Pnt(0,0,0), Vec(-1,0,0))\n",
    "right = Plane(Pnt(1,0,0), Vec(1,0,0))\n",
    "finitecyl = cyl * left * right\n",
    "geo.AddSurface(cyl, finitecyl)\n",
    "geo.NameEdge(cyl,left, \"left\")\n",
    "geo.NameEdge(cyl,right, \"right\")\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.2))\n",
    "mesh.Curve(order)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Lagrangian elements for displacement $u \\in [H^1(S)]^3$ and the rotation $\\beta \\in [H^1(S)]^3$. It might lock for small thickness $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes1 = VectorH1(mesh, order=order, dirichlet_bbnd=\"left\")\n",
    "fes = fes1*fes1\n",
    "u,beta = fes.TrialFunction()\n",
    "\n",
    "nsurf = specialcf.normal(3)\n",
    "\n",
    "thickness = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Membrane energy\n",
    "$$\n",
    "t\\|E_{tt}(u)\\|^2_{L^2(S)} \n",
    "$$\n",
    "Shear energy\n",
    "$$\n",
    "t\\int_S | \\nabla u^\\top n - \\beta |^2 \n",
    "$$\n",
    "Bending energy\n",
    "$$\n",
    "\\frac{t^3}{2}\\|\\boldsymbol{\\varepsilon}(\\beta)-\\text{Sym}(\\nabla u^\\top\\nabla\\nu)\\|^2_{L^2(S)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ptau = Id(3) - OuterProduct(nsurf,nsurf)\n",
    "Ftau = Grad(u).Trace() + Ptau\n",
    "Ctautau = Ftau.trans * Ftau\n",
    "Etautau = 0.5*(Ctautau - Ptau)\n",
    "\n",
    "eps_beta = Sym(Ptau*Grad(beta).Trace())\n",
    "gradu = Grad(u).Trace()\n",
    "ngradu = gradu.trans*nsurf\n",
    "#Average normal vector for affine geometry\n",
    "if order == 1:\n",
    "    gfn = GridFunction(fes1)\n",
    "    gfn.Set(nsurf,definedon=mesh.Boundaries(\".*\"))\n",
    "else:\n",
    "    gfn = nsurf\n",
    "\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "#membrane energy\n",
    "a += Variation( thickness*InnerProduct(Etautau, Etautau)*ds )\n",
    "#bending energy\n",
    "a += Variation( 0.5*thickness**3*InnerProduct(eps_beta-Sym(gradu.trans*Grad(gfn)),eps_beta-Sym(gradu.trans*Grad(gfn)))*ds )\n",
    "#shearing energy\n",
    "a += Variation( thickness*(ngradu-beta)*(ngradu-beta)*ds )\n",
    "\n",
    "# external force\n",
    "factor = Parameter(0.0)\n",
    "a += Variation( -thickness*factor*y*u[1]*ds )\n",
    "\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Increase the load step-wise, solve the non-linear problem by Newton's method. First and second order derivatives are computed by automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for loadstep in range(6):\n",
    "        print(\"loadstep \", loadstep)\n",
    "        factor.Set (1.5*(loadstep+1))\n",
    "        solvers.NewtonMinimization(a, gfu, printing=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(gfu.components[1], mesh, \"rotations\", deformation=gfu.components[0])\n",
    "Draw(gfu.components[0], mesh, \"disp\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Koiter shell model\n",
    "We present the method described in [<a href=\"https://www.sciencedirect.com/science/article/abs/pii/S0045794919304833\">Neunteufel and Schöberl. The Hellan-Herrmann-Johnson method for nonlinear shells. <i>Computers \\& Structures </i>, 225\n",
    "  (2019), 106109</a>]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.meshes import MakeStructuredSurfaceMesh\n",
    "thickness = 0.1\n",
    "L = 12\n",
    "W = 1\n",
    "E, nu = 1.2e6, 0\n",
    "moment = IfPos(x-L+1e-6, 1, 0)*50*pi/3\n",
    "\n",
    "mapping = lambda x,y,z : (L*x, W*y,0)\n",
    "mesh = MakeStructuredSurfaceMesh(quads=False, nx=10, ny=1, mapping=mapping)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid membrane locking Regge interpolation as in [<a href=\"https://arxiv.org/abs/1907.06232\">Neunteufel and Schöberl. Avoiding Membrane Locking with Regge Interpolation</a>] can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# False -> membrane locking\n",
    "regge = True\n",
    "order = 2\n",
    "\n",
    "fes1 = HDivDivSurface(mesh, order=order-1, discontinuous=True)\n",
    "fes2 = VectorH1(mesh, order=order, dirichletx_bbnd=\"left\", dirichlety_bbnd=\"left|bottom\", dirichletz_bbnd=\"left\")\n",
    "fes3 = HDivSurface(mesh, order=order-1, orderinner=0, dirichlet_bbnd=\"left\")\n",
    "if regge: \n",
    "    fes4 = HCurlCurl(mesh, order=order-1, discontinuous=True)\n",
    "    fes  = fes2*fes1*fes3*fes4*fes4\n",
    "    u,sigma,hyb,C,R = fes.TrialFunction()\n",
    "    sigma, hyb, C, R = sigma.Trace(), hyb.Trace(), C.Trace(), R.Operator(\"dualbnd\")\n",
    "else:\n",
    "    fes  = fes2*fes1*fes3\n",
    "    u,sigma,hyb = fes.TrialFunction()\n",
    "    sigma, hyb = sigma.Trace(), hyb.Trace()\n",
    "\n",
    "fesVF = VectorFacetSurface(mesh, order=order)\n",
    "        \n",
    "gfclamped = GridFunction(FacetSurface(mesh,order=0))\n",
    "gfclamped.Set(1,definedon=mesh.BBoundaries(\"left\"))\n",
    "\n",
    "solution = GridFunction(fes, name=\"solution\")\n",
    "averednv = GridFunction(fesVF)\n",
    "averednv_start = GridFunction(fesVF)\n",
    "        \n",
    "\n",
    "nsurf = specialcf.normal(mesh.dim)\n",
    "t     = specialcf.tangential(mesh.dim)\n",
    "nel   = Cross(nsurf, t)\n",
    "    \n",
    "Ptau    = Id(mesh.dim) - OuterProduct(nsurf,nsurf)\n",
    "Ftau    = Grad(u).Trace() + Ptau\n",
    "Ctau    = Ftau.trans*Ftau\n",
    "Etautau = 0.5*(Ctau - Ptau)\n",
    "\n",
    "nphys   = Normalize(Cof(Ftau)*nsurf)\n",
    "tphys   = Normalize(Ftau*t)\n",
    "nelphys = Cross(nphys,tphys)\n",
    "\n",
    "Hn = (u.Operator(\"hesseboundary\").trans * nphys).Reshape((3, 3))\n",
    "\n",
    "cfnphys = Normalize(Cof(Ptau+Grad(solution.components[0]))*nsurf)\n",
    "\n",
    "cfn  = Normalize(CoefficientFunction( averednv.components ))\n",
    "cfnR = Normalize(CoefficientFunction( averednv_start.components ))\n",
    "pnaverage = Normalize( cfn - (tphys*cfn)*tphys )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\sum_{T\\in \\mathcal{T}_h}\\int_{\\partial T} b\\cdot\\delta b\\,ds = \\sum_{T\\in \\mathcal{T}_h}\\int_{\\partial T} \\nu^n\\cdot\\delta b\\,ds,\\qquad \\forall \\delta b$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "averednv.Set(nsurf, dual=True, definedon=mesh.Boundaries(\".*\"))\n",
    "averednv_start.Set(nsurf, dual=True, definedon=mesh.Boundaries(\".*\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gradn = specialcf.Weingarten(3) #grad(nsurf)\n",
    "\n",
    "def MaterialNorm(mat, E, nu):\n",
    "    return E/(1-nu**2)*((1-nu)*InnerProduct(mat,mat)+nu*Trace(mat)**2)\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,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfA = BilinearForm(fes, symmetric=True, condense=True)\n",
    "bfA += Variation( (-6/thickness**3*MaterialNormInv(sigma, E, nu) \\\n",
    "                   + InnerProduct(sigma, Hn + (1-nphys*nsurf)*gradn))*ds ).Compile()\n",
    "if regge:\n",
    "    bfA += Variation( 0.5*thickness*MaterialNorm(C, E, nu)*ds )\n",
    "    bfA += Variation( InnerProduct(C-Etautau, R)*ds(element_vb=BND) )\n",
    "    bfA += Variation( InnerProduct(C-Etautau, R)*ds(element_vb=VOL) )\n",
    "else:\n",
    "    bfA += Variation( 0.5*thickness*MaterialNorm(Etautau, E, nu)*ds )\n",
    "bfA += Variation( -(acos(nel*cfnR)-acos(nelphys*pnaverage)-hyb*nel)*(sigma*nel)*nel*ds(element_boundary=True) ).Compile()\n",
    "par = Parameter(0.0)\n",
    "bfA += Variation( -par*moment*(hyb*nel)*ds(element_boundary=True) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par.Set(0.1)\n",
    "averednv.Set((1 - gfclamped) * cfnphys + gfclamped * nsurf,\n",
    "            definedon=mesh.Boundaries(\".*\"),\n",
    "            dual=True,\n",
    "        )\n",
    "\n",
    "with TaskManager():\n",
    "    solvers.Newton(bfA, solution, inverse=\"sparsecholesky\", maxerr=1e-10, maxit=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(solution.components[0], mesh, \"disp\", deformation=solution.components[0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numsteps=10\n",
    "with TaskManager():\n",
    "    for steps in range(1,numsteps):\n",
    "        par.Set((steps+1)/numsteps)\n",
    "        print(\"Loadstep =\", steps+1, \", F/Fmax =\", (steps+1)/numsteps*100, \"%\")\n",
    "        \n",
    "        averednv.Set((1 - gfclamped) * cfnphys + gfclamped * nsurf,\n",
    "            definedon=mesh.Boundaries(\".*\"),\n",
    "            dual=True,\n",
    "        )\n",
    "        \n",
    "        (res,numit) = solvers.Newton(bfA, solution, inverse=\"sparsecholesky\", printing=False, maxerr=2e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(solution.components[0], mesh, \"disp\", deformation=solution.components[0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}