NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dielectric shell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "import numpy as np\n",
    "from netgen.csg import *\n",
    "\n",
    "from newtonmethod import NewtonWithLinesearch\n",
    "\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "# SetNumThreads(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thickness = 0.5\n",
    "RAD = 20\n",
    "mu = 1e-2\n",
    "\n",
    "\n",
    "Voltage = 5000\n",
    "epsilon0 = 8.854e-12\n",
    "epsilon = 5*epsilon0\n",
    "N_chain = 5\n",
    "capacity = epsilon/(thickness)\n",
    "charge = Voltage/(thickness)*capacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 2\n",
    "\n",
    "\n",
    "BCclamp = \"sym.|bottom\"\n",
    "geo          = CSGeometry()\n",
    "sphere       = Sphere(Pnt(0,0,0), RAD)\n",
    "bot          = Plane(Pnt(0,0,0), Vec(0,0,-1)).maxh(RAD*np.sin(np.pi/2/50))\n",
    "bot1          = Plane(Pnt(0,0,RAD*np.sin(np.pi/2/60)), Vec(0,0,-1)).maxh(RAD*np.sin(np.pi/2/50))\n",
    "bot2          = Plane(Pnt(0,0,RAD*np.sin(np.pi/2/30)), Vec(0,0,-1)).maxh(RAD*np.sin(np.pi/2/50))\n",
    "bot3          = Plane(Pnt(0,0,RAD*np.sin(np.pi/2/18)), Vec(0,0,-1)).maxh(RAD*np.sin(np.pi/2/50))\n",
    "bot4          = Plane(Pnt(0,0,RAD*np.sin(np.pi/2/12)), Vec(0,0,-1)).maxh(RAD*np.sin(np.pi/2/50))\n",
    "# top          = Plane(Pnt(0,0,W), Vec(0,0,1))\n",
    "planex = Plane(Pnt(0,0,0), Vec(-1,0,0))\n",
    "planey = Plane(Pnt(0,0,0), Vec(0,-1,0))\n",
    "finitesphere1 = sphere * bot * planex * planey - bot1\n",
    "finitesphere2 = sphere * bot1 * planex * planey - bot2\n",
    "finitesphere3 = sphere * bot2 * planex * planey - bot3\n",
    "finitesphere4 = sphere * bot3 * planex * planey - bot4\n",
    "finitesphere5 = sphere * bot4 * planex * planey\n",
    "\n",
    "\n",
    "geo.AddSurface(sphere, finitesphere1.bc(\"surface\"))\n",
    "geo.AddSurface(sphere, finitesphere2.bc(\"surface\"))\n",
    "geo.AddSurface(sphere, finitesphere3.bc(\"surface\"))\n",
    "geo.AddSurface(sphere, finitesphere4.bc(\"surface\"))\n",
    "geo.AddSurface(sphere, finitesphere5.bc(\"surface\"))\n",
    "geo.NameEdge(sphere,bot, \"bottom\")\n",
    "geo.NameEdge(sphere,planex, \"symx\")\n",
    "geo.NameEdge(sphere,planey, \"symy\")\n",
    "\n",
    "geo.AddPoint(Pnt(RAD,0,0), \"pntload_x\")\n",
    "geo.AddPoint(Pnt(0,RAD,0), \"pntload_y\")\n",
    "\n",
    "\n",
    "from netgen.meshing import MeshingParameters\n",
    "maxhrel = 9\n",
    "mp = MeshingParameters(maxh=RAD/maxhrel, grading=0.1)\n",
    "mesh = Mesh(geo.GenerateMesh(mp))\n",
    "\n",
    "mesh.Curve(order)\n",
    "\n",
    "Draw(mesh)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_mom = HDivDivSurface(mesh, order=order-1, discontinuous=True)\n",
    "fes_u = VectorH1(mesh, order=order, dirichletx_bbnd=\"symx|bottom\", dirichlety_bbnd=\"symy|bottom\", dirichletz_bbnd=\"bottom\")\n",
    "fes_hyb = HDivSurface(mesh, order=order-1, orderinner=0, dirichlet_bbnd=BCclamp)\n",
    "fes_curlcurl = HCurlCurl(mesh, order=order-1, discontinuous=True, print=True)\n",
    "fes_l2 = SurfaceL2(mesh, order=order-1)\n",
    "\n",
    "fes  = fes_u*fes_hyb*fes_curlcurl*fes_curlcurl*fes_curlcurl*fes_mom*fes_l2**2 * fes_l2**2 * fes_l2\n",
    "\n",
    "#%%\n",
    "par = Parameter(0.0)\n",
    "\n",
    "solution = GridFunction(fes, name=\"solution\")\n",
    "solution_0 = GridFunction(fes)\n",
    "u = solution.components[0]\n",
    "eps = solution.components[2]\n",
    "kappa = solution.components[4]\n",
    "mom = solution.components[5]\n",
    "\n",
    "\n",
    "Nsurf = specialcf.normal(mesh.dim) # surface normal N\n",
    "t     = specialcf.tangential(mesh.dim)\n",
    "nel   = Cross(Nsurf, t) # in-plane edge normal\n",
    "\n",
    "gradN = specialcf.Weingarten(3) \n",
    "\n",
    "\n",
    "A = Id(mesh.dim) - OuterProduct(Nsurf,Nsurf)\n",
    "\n",
    "cfnphys = Normalize(Cof(A+Grad(u))*Nsurf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fesVF = VectorFacetSurface(mesh, order=order)\n",
    "averednv = GridFunction(fesVF)\n",
    "averednv_start = GridFunction(fesVF)\n",
    "\n",
    "\n",
    "fesF = FacetSurface(mesh, order=0)\n",
    "gfclamped = GridFunction(fesF)\n",
    "gfclamped.Set(1,definedon=mesh.BBoundaries(BCclamp))\n",
    "\n",
    "n_ = fesVF.TrialFunction()\n",
    "n_.Reshape((3,))\n",
    "bfF = BilinearForm(fesVF, symmetric=True)\n",
    "bfF += Variation( (0.5*n_*n_ - ((1-gfclamped)*cfnphys+gfclamped*Nsurf)*n_)*ds(element_boundary=True))\n",
    "\n",
    "\n",
    "def ComputeAveredNV(averednv):\n",
    "    rf = averednv.vec.CreateVector()\n",
    "    bfF.Apply(averednv.vec, rf)\n",
    "    bfF.AssembleLinearization(averednv.vec)\n",
    "    invF = bfF.mat.Inverse(fesVF.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "    averednv.vec.data -= invF*rf\n",
    "\n",
    "\n",
    "ComputeAveredNV(averednv)\n",
    "ComputeAveredNV(averednv_start)\n",
    "\n",
    "cfn  = Normalize(CoefficientFunction( averednv.components ))\n",
    "cfnR = Normalize(CoefficientFunction( averednv_start.components ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B0 = GridFunction(fes_curlcurl)\n",
    "def ComputeB(B0):\n",
    "    # for precomputing B0\n",
    "    fes3B0 = HDivSurface(mesh, order=order-1, orderinner=0, dirichlet_bbnd=\".*\")\n",
    "    fes_B0  = fes3B0*fes_curlcurl*fes_mom\n",
    "    sol = GridFunction(fes_B0)\n",
    "\n",
    "    hyb_, B_, mom_ = fes_B0.TrialFunction()\n",
    "    hyb_, B_, mom_ = hyb_.Trace(), B_.Trace(), mom_.Trace()\n",
    "\n",
    "    bfA0 = BilinearForm(fes_B0, symmetric=True)\n",
    "    bfA0 += Variation( InnerProduct(B_, B_)*ds).Compile()\n",
    "    bfA0 += Variation( (-InnerProduct(mom_, B_ + gradN))*ds ).Compile()\n",
    "    bfA0 += Variation( (acos(nel*cfnR)-np.pi/2-hyb_*nel)*(mom_*nel)*nel*ds(element_boundary=True ) ).Compile()\n",
    "\n",
    "    sol.vec[:] = 0\n",
    "    NewtonWithLinesearch(bfA0, sol.vec, inverse=\"sparsecholesky\")\n",
    "        \n",
    "    B0.vec.data = sol.components[1].vec\n",
    "\n",
    "ComputeB(B0)\n",
    "Draw(B0, mesh, \"B\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_, hyb_, eps_, R_, kappa_, mom_ = fes.TrialFunction()[:6]\n",
    "mom_, hyb_, eps_, R_, kappa_ = mom_.Trace(), hyb_.Trace(), eps_.Trace(), R_.Operator(\"dualbnd\"), kappa_.Trace()\n",
    "lam_, pres_, E1_ = fes.TrialFunction()[6:]\n",
    "lam_, pres_, E1_ = lam_.Trace(), pres_.Trace(), E1_.Trace()\n",
    "\n",
    "Fsurf_    = grad(u_).Trace() + A\n",
    "Csurf_    = Fsurf_.trans*Fsurf_\n",
    "epssurf_ = 0.5*(Csurf_ - A)\n",
    "\n",
    "nphys   = Normalize(Cof(Fsurf_)*Nsurf) #  normal of deformed surface\n",
    "tphys   = Normalize(Fsurf_*t)\n",
    "nelphys = Cross(nphys,tphys) # in-plane edge normal of deformed surface\n",
    "\n",
    "Hn_ = CoefficientFunction( (u_.Operator(\"hesseboundary\").trans*nphys), dims=(3,3) )\n",
    "pnaverage = Normalize( cfn - (tphys*cfn)*tphys )\n",
    "\n",
    "\n",
    "bfA = BilinearForm(fes, symmetric=True, condense=True)\n",
    "gausspoints = [(-np.sqrt(3/5), 5/9), (0,8/9), ( np.sqrt(3/5), 5/9)]\n",
    "\n",
    "for (zi, wi) in gausspoints:\n",
    "    zeta = thickness/2*zi\n",
    "    weightdet = wi*thickness/2\n",
    "    lamz = (lam_[0]+1)*zeta + lam_[1]/2*zeta*zeta\n",
    "    lamzp = lam_[0]+1 + lam_[1]*zeta\n",
    "    presz = pres_[0] + pres_[1]*zeta\n",
    "\n",
    "    C_ = A + 2*(eps_ + lamz*kappa_ - lamz*B0 + zeta*B0 + zeta*eps_*B0 + zeta*B0*eps_)\n",
    "\n",
    "    I_C = Trace(C_)\n",
    "    CNN = lamzp**2\n",
    "    III_C2 = (Cof(C_)*Nsurf)*Nsurf\n",
    "    FB = (-zeta*B0) + A\n",
    "    III_FB = (Cof(FB)*Nsurf)*Nsurf\n",
    "\n",
    "    detC = III_C2*CNN\n",
    "    I_barC = detC**(-1/3)*(I_C + CNN)\n",
    "    \n",
    "    ArrudaBoyce = 0.5*mu*(I_barC - 3) + mu/20/N_chain*(I_barC**2 - 9) + \\\n",
    "        11*mu/1050/N_chain**2*(I_barC**3 - 3**3) + \\\n",
    "        19*mu/7000/N_chain**3*(I_barC**4 - 3**4) + \\\n",
    "        519*mu/673750/N_chain**4*(I_barC**5 - 3**5)\n",
    "    bfA += Variation(weightdet*III_FB*(ArrudaBoyce)*ds).Compile()\n",
    "    Efield_r = -par*Voltage/thickness + E1_*zeta\n",
    "    bfA += Variation(-0.5*epsilon * Efield_r**2*III_FB/CNN*weightdet*ds).Compile()\n",
    "    bfA += Variation( (sqrt(III_C2*CNN)-1) * presz *III_FB*weightdet*ds)\n",
    "\n",
    "\n",
    "bfA += Variation( (InnerProduct(mom_, kappa_ + Hn_ + (1-nphys*Nsurf)*gradN))*ds ).Compile()\n",
    "bfA += Variation( InnerProduct(eps_-epssurf_, R_)*ds(element_vb=BND) )\n",
    "bfA += Variation( InnerProduct(eps_-epssurf_, R_)*ds(element_vb=VOL) )\n",
    "bfA += Variation( -(acos(nel*cfnR)-acos(nelphys*pnaverage)-hyb_*nel)*(mom_*nel)*nel*ds(element_boundary=True) ).Compile()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ndof_coupling = sum(fes.FreeDofs(True))\n",
    "print(f\"coupling dofs = {ndof_coupling}\")\n",
    "\n",
    "\n",
    "# sceneu = Draw(u, mesh, \"disp\", deformation=u)\n",
    "sceneur = Draw(u*Nsurf, mesh, \"dispr\", deformation=u)\n",
    "\n",
    "# scenem = Draw(mom, mesh, \"m\", deformation=u)\n",
    "scenekappa = Draw(kappa, mesh, \"kappa\", deformation=u)\n",
    "\n",
    "#%%\n",
    "\n",
    "glist = np.concatenate((np.linspace(0,0.8,4,endpoint=False),\n",
    "                        np.linspace(0.8,1,4, endpoint=False),\n",
    "                        np.linspace(1.0,1.03,3, endpoint=False),\n",
    "                        np.linspace(1.03,1.05,5, endpoint=False),\n",
    "                        np.linspace(1.05,1.06,20, endpoint=False),\n",
    "                        np.linspace(1.06,1.065,11, endpoint=True)))\n",
    "\n",
    "\n",
    "\n",
    "#%%\n",
    "counter = 0\n",
    "\n",
    "lastloadpar = glist[0]\n",
    "with TaskManager():\n",
    "    for loadpar in glist[1:]:\n",
    "        solution_0.vec.data = solution.vec\n",
    "        par.Set(loadpar)\n",
    "        print(\"Loadstep value =\", par.Get())\n",
    "        print(\"Applied voltage =\", Voltage*par.Get())\n",
    "\n",
    "        # ComputeAveredNV(averednv)\n",
    "\n",
    "        (err,numit) = NewtonWithLinesearch(a=bfA, x=solution.vec, abserror=1e-7, maxnewton=40, factor=1.03)\n",
    "\n",
    "\n",
    "\n",
    "        if err > 1e-3:\n",
    "            print(\"Newton method did not converge\")\n",
    "            solution.vec.data = solution_0.vec\n",
    "            par.Set(lastloadpar)\n",
    "            break\n",
    "\n",
    "        poslist = []\n",
    "\n",
    "        for theta in np.linspace(0,np.pi/2,101):\n",
    "            ur = u*Nsurf\n",
    "            pos = ur(mesh(RAD*np.cos(np.pi/2/18)*np.cos(theta), RAD*np.cos(np.pi/2/18)*np.sin(theta), RAD*np.sin(np.pi/2/18), BND))\n",
    "            poslist.append(pos)\n",
    "            # print(pos, end=\"\\t\", file=outfile)\n",
    "\n",
    "        print(\"min/max\", np.min(poslist), np.max(poslist))\n",
    "\n",
    "        sceneur.Redraw()\n",
    "        scenekappa.Redraw()\n",
    "        lastloadpar = loadpar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}