{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f987fdcb-5bb6-49a5-b68f-917c95a43f55",
   "metadata": {},
   "source": [
    "### Dirichlet boundary conditions¶\n",
    "This example shows in depth into the mechanisms required to solve a model with non-homogenous Dirichlet BCs (see https://docu.ngsolve.org/nightly/i-tutorials/unit-1.3-dirichlet/dirichlet.html). In Detail: \n",
    "\n",
    "  * The model is a 3D bar.\n",
    "  * Combinations of homogenous (fix) and non-homogenous Dirichlet BCs (moving BC) are tested.\n",
    "  * The standard and CGSolver can be used.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2fc5ef92-754b-4e0c-a7c0-1d49e23bb0d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Left', 'Right', 'Bottom', 'Top', 'Front', 'Back']\n"
     ]
    }
   ],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "# Variables\n",
    "E = 210000\n",
    "nu = 0.3\n",
    "\n",
    "# Dirichlet Boundaries (homogenous)\n",
    "Face_Dirichlet = 'Left'\n",
    "# Displacement boundary (non-homogenous)\n",
    "FaceDisplacement = 'Right'\n",
    "\n",
    "# Degree of Element\n",
    "EleDeg = 2\n",
    "\n",
    "# Mesh the geometry and get the boundaries --------------------------------------------------\n",
    "box = Box((0,-10,-10), (200,10,10))\n",
    "box.faces.Min(X).name=\"Left\"\n",
    "box.faces.Max(X).name=\"Right\"\n",
    "box.faces.Min(Y).name=\"Bottom\"\n",
    "box.faces.Max(Y).name=\"Top\"\n",
    "box.faces.Min(Z).name=\"Front\"\n",
    "box.faces.Max(Z).name=\"Back\"\n",
    "\n",
    "geo = OCCGeometry(box)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)\n",
    "boundaries = list(mesh.GetBoundaries())\n",
    "print(boundaries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "93ead963-74ca-4c46-9070-6724e5fe5c3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bede6a47bade410fa1b8d73d0489f8fa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='400px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Draw(mesh,  height=\"400px\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fd30d120-0d1a-4ede-b19b-a697471614bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stress Function ----------------------------------------------------------------------------\n",
    "mu = E / 2 / (1 + nu)\n",
    "lam = E * nu / ((1 + nu) * (1 - 2 * nu))\n",
    "\n",
    "def Stress(strain):\n",
    "    return 2 * mu * strain + lam * Trace(strain) * Id(3)\n",
    "\n",
    "# FE Space ---------------------------------------------------------------------------------\n",
    "# fix - fix BCs\n",
    "fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichlet, dirichlety=Face_Dirichlet, dirichletz=Face_Dirichlet, dirichlet=FaceDisplacement)\n",
    "\n",
    "# fix - fix BCs\n",
    "#fes = VectorH1(mesh, order=EleDeg,  dirichlet='Left|Right')\n",
    "\n",
    "# pin - fix BCs \n",
    "#fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichlet, dirichlet=FaceDisplacement)\n",
    "\n",
    "# fix - pin BCs\n",
    "#fes = VectorH1(mesh, order=EleDeg,  dirichlet=Face_Dirichlet, dirichletx=FaceDisplacement)\n",
    "\n",
    "# pin - pin BCs\n",
    "#fes = VectorH1(mesh, order=EleDeg,  dirichletx='Left|Right', dirichlety='Top', dirichletz='Back')\n",
    "\n",
    "u, v = fes.TnT()\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "92ff6a6a-db79-47de-ad6c-9e9a517eeedd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a6d449f007b4f179ce3856ec5ce335c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.26…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Apply displacement boundary condition -----------------------------------------------------\n",
    "displacement = CF((1, 0, 0))\n",
    "gfu = GridFunction(fes)\n",
    "gfu.Set(displacement, definedon=mesh.Boundaries(FaceDisplacement))\n",
    "Draw(gfu);\n",
    "\n",
    "# Assemble a ----------------------------------------------------------------------------\n",
    "with TaskManager():\n",
    "    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile() * dx)\n",
    "    #pre = Preconditioner(a, \"bddc\")  # for CGSolver\n",
    "    a.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6ec71cb5-5ed5-45bb-b74e-fb6642dc6fc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Non-zero force needed -> e.g. gravity along y \n",
    "force = CF((0, 0.0001, 0))\n",
    "f = LinearForm(force*v*dx).Assemble()\n",
    "r = f.vec - a.mat * gfu.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "17fea4db-dd6e-45e4-9b3a-fca04903c8c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d7a4a88011d944adb178f297f3e38606",
       "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"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run the Solver --------------------------------------------------------------------------\n",
    "\n",
    "# CGSolver\n",
    "#from ngsolve.krylovspace import CGSolver\n",
    "#inv = CGSolver(a.mat, pre, printrates='\\r', tol=1e-6, maxiter=500)\n",
    "\n",
    "# standard solver\n",
    "inv = a.mat.Inverse(freedofs=fes.FreeDofs())\n",
    "\n",
    "# solve \n",
    "gfu.vec.data += inv * r\n",
    "\n",
    "# draw displacements\n",
    "Draw (gfu, mesh, height=\"500px\", deformation=True, scale=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2cafbf12-d75e-4e97-b4fd-376e00b6878f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b8b8fa28729d4c56ae2fb14aa0e625ae",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='400px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# finite element stress\n",
    "fesStress = MatrixValued(H1(mesh, order=EleDeg), symmetric=True)\n",
    "# map to grid funtion \n",
    "gfStress = GridFunction(fesStress)\n",
    "gfStress.Interpolate(Stress(Sym(Grad(gfu))))\n",
    "Draw (Norm(gfStress),mesh, height=\"400px\", deformation=1*gfu, draw_vol=False)"
   ]
  }
 ],
 "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.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}