Nonlinear_elasticity_gels_MASU.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Free swelling using geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from ngsolve import *\n",
    "from netgen.csg import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd9b3f685af94f94a16524cd8dcdeb56",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('left', 'bot', 'default', 'default', 'back', 'default')\n"
     ]
    }
   ],
   "source": [
    "# Geometry of 1/8 \n",
    "L = 90\n",
    "d = 23.5\n",
    "L3 = 3.0\n",
    "\n",
    "left  = Plane (Pnt(0,0,0), Vec(-1,0,0) ).bc('left')\n",
    "right = Plane (Pnt(0.5*L,0,0), Vec( 1,0,0) )\n",
    "bot = Plane (Pnt(0,0,0), Vec(0,-1,0) ).bc('bot')\n",
    "top  = Plane (Pnt(0,0.5*d,0), Vec(0, 1,0) )\n",
    "back   = Plane (Pnt(0,0,0), Vec(0,0,-1) ).bc('back')\n",
    "front   = Plane (Pnt(0,0,0.5*L3), Vec(0,0, 1) )\n",
    "\n",
    "brick = left * right * front * back * bot * top\n",
    "\n",
    "geo = CSGeometry()\n",
    "geo.Add (brick)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=1))\n",
    "Draw(mesh)\n",
    "print(mesh.GetBoundaries())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma target: 0.001145522071527073\n",
      "shear modulus target: 0.15647831497059816\n",
      "lo que tiene que crecer: 11.515\n",
      "desplazamineto esperados: [22.05    5.7575  0.735 ]\n"
     ]
    }
   ],
   "source": [
    "global phi0\n",
    "#phi0 = 0.2035 # or 0.3\n",
    "phi0 = 0.3\n",
    "\n",
    "global chi \n",
    "chi = 0.45\n",
    "\n",
    "global entropic_unit\n",
    "entropic_unit=136.6   #This is (k_B * T)/(V_m), measured in Megapascals\n",
    "\n",
    "global G    #It is the shear modulus. Will be defined in the next script (next jupy cell).\n",
    "\n",
    "phi = lambda J: phi0/J\n",
    "dH = lambda J: (1-1/1000) * phi(J) + np.log(1-phi(J)) + chi * phi(J)**2\n",
    "\n",
    "#gammafun = lambda lamb: -dH(lamb)/lamb # bonded\n",
    "gammafun = lambda lamb: -dH(lamb**3)*lamb\n",
    "\n",
    "lambda_target = 1.49\n",
    "gamma_target = gammafun(lambda_target)\n",
    "G_target = gamma_target*entropic_unit\n",
    "print(\"gamma target:\", gamma_target)\n",
    "print(\"shear modulus target:\", G_target)\n",
    "\n",
    "print(\"lo que tiene que crecer:\", d*(lambda_target-1))\n",
    "G=G_target\n",
    "print( \"desplazamineto esperados:\", (lambda_target-1)/2*np.array([L,d,L3]) )\n",
    "#v = np.linspace(1,2,100)\n",
    "#plt.plot(v,gammafun(v)*entropic_unit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.001145522071527073\n"
     ]
    }
   ],
   "source": [
    "def F(u):\n",
    "    return Id(mesh.dim) + Grad(u)\n",
    "\n",
    "global gamma\n",
    "gamma = G/entropic_unit\n",
    "print(gamma)\n",
    "\n",
    "def Gels (F):\n",
    "    J = Det(F)\n",
    "    phi = phi0/J\n",
    "    H = (J - phi0)*log(1-phi)  + phi0 * chi*(1-phi) + phi0/1000*log(phi)\n",
    "    C = F.trans * F\n",
    "    return 0.5*gamma*(Trace(C)) + H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finite element space with slipp boundary conditions on back|bot|left\n",
    "V = VectorH1(mesh, order=2, dirichletx = 'left', dirichlety='bot', dirichletz='back')\n",
    "\n",
    "# Construction of bilinear form\n",
    "u  = V.TrialFunction()\n",
    "I = Id(mesh.dim)\n",
    "F = I + Grad(u)\n",
    "a = BilinearForm(V, symmetric=True)\n",
    "\n",
    "# hydrogel model\n",
    "a += Variation(  Gels (F).Compile() * dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveNonlinearMinProblem(a,gfu,tol=1e-08,maxits=250, scenes=None):\n",
    "    res = gfu.vec.CreateVector()\n",
    "    du  = gfu.vec.CreateVector()\n",
    "\n",
    "    for it in range(maxits):\n",
    "        print (\"Newton iteration {:3}\".format(it),end=\", \")\n",
    "        print (\"energy = {:16}\".format(a.Energy(gfu.vec)),end=\"\")\n",
    "        #solve linearized problem:\n",
    "        a.Apply (gfu.vec, res)\n",
    "        a.AssembleLinearization (gfu.vec)\n",
    "        inv = a.mat.Inverse(V.FreeDofs())\n",
    "        alpha = 5e-1\n",
    "        du.data = alpha * inv * res\n",
    "\n",
    "        #update iteration\n",
    "        gfu.vec.data -= du\n",
    "\n",
    "        #stopping criteria\n",
    "        stopcritval = sqrt(abs(InnerProduct(du,res)))\n",
    "        print (\"<A u\",it,\", A u\",it,\">_{-1}^0.5 = \", stopcritval)\n",
    "        if stopcritval < tol:\n",
    "            break\n",
    "        for sc in scenes:\n",
    "            sc.Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bdcebefb3b1b4606ba65c154882822a4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton iteration   0, energy = -121.99481492764417<A u 0 , A u 0 >_{-1}^0.5 =  1.581255529712478\n",
      "Newton iteration   1, energy = -123.97351994566714<A u 1 , A u 1 >_{-1}^0.5 =  1.2859745704927337\n",
      "Newton iteration   2, energy = -125.27881870605432<A u 2 , A u 2 >_{-1}^0.5 =  1.0204961804833017\n",
      "Newton iteration   3, energy = -126.09747529045949<A u 3 , A u 3 >_{-1}^0.5 =  0.7795333522749901\n",
      "Newton iteration   4, energy = -126.57231306714256<A u 4 , A u 4 >_{-1}^0.5 =  0.5633469922079162\n",
      "Newton iteration   5, energy = -126.8183495534079<A u 5 , A u 5 >_{-1}^0.5 =  0.37888040355189156\n",
      "Newton iteration   6, energy = -126.92865455723854<A u 6 , A u 6 >_{-1}^0.5 =  0.23527236168219756\n",
      "Newton iteration   7, energy = -126.970837018616<A u 7 , A u 7 >_{-1}^0.5 =  0.13567495759027837\n",
      "Newton iteration   8, energy = -126.98477507277244<A u 8 , A u 8 >_{-1}^0.5 =  0.07388725225427868\n",
      "Newton iteration   9, energy = -126.98889132710389<A u 9 , A u 9 >_{-1}^0.5 =  0.038743318375331084\n",
      "Newton iteration  10, energy = -126.99002027846409<A u 10 , A u 10 >_{-1}^0.5 =  0.019866986796776526\n",
      "Newton iteration  11, energy = -126.99031673062437<A u 11 , A u 11 >_{-1}^0.5 =  0.010063782826685226\n",
      "Newton iteration  12, energy = -126.99039274636087<A u 12 , A u 12 >_{-1}^0.5 =  0.005065327800197487\n",
      "Newton iteration  13, energy = -126.99041199666252<A u 13 , A u 13 >_{-1}^0.5 =  0.0025411348468145794\n",
      "Newton iteration  14, energy = -126.99041684058925<A u 14 , A u 14 >_{-1}^0.5 =  0.0012726993917459744\n",
      "Newton iteration  15, energy = -126.99041805552567<A u 15 , A u 15 >_{-1}^0.5 =  0.0006368844828935932\n",
      "Newton iteration  16, energy = -126.99041835975602<A u 16 , A u 16 >_{-1}^0.5 =  0.0003185761635775158\n",
      "Newton iteration  17, energy = -126.99041843587608<A u 17 , A u 17 >_{-1}^0.5 =  0.0001593215905571567\n",
      "Newton iteration  18, energy = -126.99041845491377<A u 18 , A u 18 >_{-1}^0.5 =  7.966917600506443e-05\n",
      "Newton iteration  19, energy = -126.9904184596743<A u 19 , A u 19 >_{-1}^0.5 =  3.983668362428844e-05\n",
      "Newton iteration  20, energy = -126.99041846086448<A u 20 , A u 20 >_{-1}^0.5 =  1.9918865770881715e-05\n",
      "Newton iteration  21, energy = -126.99041846116195<A u 21 , A u 21 >_{-1}^0.5 =  9.959563886678913e-06\n",
      "Newton iteration  22, energy = -126.99041846123652<A u 22 , A u 22 >_{-1}^0.5 =  4.9798146912265985e-06\n",
      "Newton iteration  23, energy = -126.99041846125503<A u 23 , A u 23 >_{-1}^0.5 =  2.4899155334282167e-06\n",
      "Newton iteration  24, energy = -126.99041846125965<A u 24 , A u 24 >_{-1}^0.5 =  1.2449598128795222e-06\n",
      "Newton iteration  25, energy = -126.99041846126083<A u 25 , A u 25 >_{-1}^0.5 =  6.224804183813055e-07\n",
      "Newton iteration  26, energy = -126.99041846126111<A u 26 , A u 26 >_{-1}^0.5 =  3.112403356523248e-07\n",
      "Newton iteration  27, energy = -126.99041846126117<A u 27 , A u 27 >_{-1}^0.5 =  1.5562020312226921e-07\n",
      "Newton iteration  28, energy = -126.99041846126119<A u 28 , A u 28 >_{-1}^0.5 =  7.781010787623135e-08\n",
      "Newton iteration  29, energy = -126.99041846126116<A u 29 , A u 29 >_{-1}^0.5 =  3.890505690861817e-08\n",
      "Newton iteration  30, energy = -126.9904184612612<A u 30 , A u 30 >_{-1}^0.5 =  1.945252858667614e-08\n",
      "Newton iteration  31, energy = -126.99041846126119<A u 31 , A u 31 >_{-1}^0.5 =  9.726266036205506e-09\n"
     ]
    }
   ],
   "source": [
    "# Solve\n",
    "\n",
    "gfu = GridFunction(V)\n",
    "gfu.vec[:] = 0\n",
    "\n",
    "scene = Draw (gfu, mesh, \"u\", deformation=True)\n",
    "\n",
    "# scene0 = Draw (gfu.components[0], mesh, \"u\", deformation=True)\n",
    "# scene1 = Draw (gfu.components[1], mesh, \"u\", deformation=True)\n",
    "\n",
    "# scene2 = Draw (gfu.components[2], mesh, \"u\", deformation=True)\n",
    "\n",
    "SetVisualization (deformation=True)\n",
    "\n",
    "res = gfu.vec.CreateVector()\n",
    "du = gfu.vec.CreateVector()\n",
    "\n",
    "# scenes = [scene, scene0, scene1, scene2]\n",
    "scenes = [scene]\n",
    "SolveNonlinearMinProblem(a,gfu, scenes=scenes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'result_free_swelling_sub3'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VTKOutput object\n",
    "vtk = VTKOutput(ma=mesh,\n",
    "                coefs=[gfu],\n",
    "                names = [\"displacement\"],\n",
    "                filename=\"result_free_swelling_sub3\",\n",
    "                subdivision=3)\n",
    "# Exporting the results:\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'RUN_Gel_free_swelling'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[9], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mRUN_Gel_free_swelling\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m gel_free_swelling, Solve_gel_free_swelling\n\u001b[0;32m      2\u001b[0m GEL \u001b[38;5;241m=\u001b[39m gel_free_swelling()\n\u001b[0;32m      3\u001b[0m GEL\u001b[38;5;241m.\u001b[39mgeometry()\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'RUN_Gel_free_swelling'"
     ]
    }
   ],
   "source": [
    "from RUN_Gel_free_swelling import gel_free_swelling, Solve_gel_free_swelling\n",
    "GEL = gel_free_swelling()\n",
    "GEL.geometry()\n",
    "\n",
    "modelling = Solve_gel_free_swelling(GEL)\n",
    "modelling.Solve()\n",
    "Draw(modelling.gfu)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "interpreter": {
   "hash": "e51f2c4ff3136256e33bf4519c50b38d40865c48a21fffea67573bae2ea1b8ed"
  },
  "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.10.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}