{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## NGSolve Model Templates\n",
    "* A ready to use collection of physical models based on NGSolve\n",
    "* Solid mechanics, Fluid Dynamics, Electromagnetics ...\n",
    "* Easy to combine models to solve multi-physics problems\n",
    "* Useful as basis for optimization, model order reduction, UQ, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "... work in progress"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2001 - 2015: C++ code only \n",
    "* A lot of special files and #ifdef SABINE\n",
    "* Exchange C++ files and snippets, had to reconfigure and recompile \n",
    "\n",
    "## 2015 - xxxx: Python frontend\n",
    "* Unified C++ code (with some extensions)\n",
    "* Python formulation much shorter and expressive, good i-Tutorials\n",
    "* Exchange Python files and snippets \n",
    "* Still a lot of expert knowledge required for advanced methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "import whatever we need "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import math, time\n",
    "\n",
    "from ngsolve import *\n",
    "import netgen.geom2d\n",
    "import netgen.meshing\n",
    "import netgen.gui\n",
    "\n",
    "from ngsolve.internal import visoptions\n",
    "\n",
    "from ngs_templates.Elasticity import * \n",
    "from ngs_templates.NavierStokes import *\n",
    "from ngs_templates.Transport import *\n",
    "from ngs_templates.ConvectionDiffusion import * \n",
    "\n",
    "print (\"gui ist up\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Elasticity Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "geo = netgen.geom2d.SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (1, 0.1), bcs = (\"bottom\", \"right\", \"top\", \"left\"))\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "\n",
    "model = Elasticity(mesh=mesh, materiallaw=HookeMaterial(E=200,nu=0.2),\n",
    "                    dirichlet=\"left\",\n",
    "                    boundaryforce = { \"right\" : (0,0.2) },\n",
    "                    nonlinear=True, order=4)\n",
    "\n",
    "ngsglobals.msg_level = 0\n",
    "model.Solve()\n",
    "\n",
    "Draw (model.displacement, reset=True)\n",
    "Draw (model.stress[0,0], mesh, \"stress_xx\")\n",
    "SetVisualization (deformation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "dispright = Parameter(0)\n",
    "model = Elasticity(mesh=mesh, materiallaw=HookeMaterial(E=200,nu=0.2), \\\n",
    "            boundarydisplacement = { \"left\" : (0,0), \"right\" : (0,dispright)},\n",
    "            nonlinear=True, order=4)\n",
    "\n",
    "\n",
    "Draw (model.displacement, reset=True) \n",
    "Draw (model.stress[0,0], mesh, \"stress_xx\")\n",
    "SetVisualization (deformation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dispright.Set (dispright.Get()+0.1)\n",
    " \n",
    "ngsglobals.msg_level = 0\n",
    "model.Solve()\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# mesh-file available from:\n",
    "# https://nemesis.asc.tuwien.ac.at/index.php/s/fDdc6gn5bRRScJK\n",
    "url = 'https://nemesis.asc.tuwien.ac.at/index.php/s/fDdc6gn5bRRScJK/download?path=%2F&files=frame.unv&downloadStartSecret=76ip0ltng4u'\n",
    "import os, requests\n",
    "if not os.path.exists('frame.unv'):\n",
    "    import requests\n",
    "    r = requests.get(url, allow_redirects=True)\n",
    "    open('frame.unv', 'wb').write(r.content)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh(netgen.meshing.ImportMesh(\"frame.unv\"))\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model = Elasticity(mesh=mesh, materiallaw=HookeMaterial(200,0.2), \\\n",
    "                       dirichlet=\"holes\",\n",
    "                       boundaryforce = { \"right\" : (0,0,1), \"left\" : (0,0,-1) },\n",
    "                       nonlinear=False, order=2)\n",
    "\n",
    "with TaskManager():\n",
    "    model.Solve()\n",
    "\n",
    "Draw (model.displacement, mesh, \"disp\", reset=True)\n",
    "Draw (model.stress, mesh, \"stress\")\n",
    "\n",
    "SetVisualization (deformation=True)\n",
    "\n",
    "myfes = H1(mesh, order=2)\n",
    "normstress = GridFunction(myfes)\n",
    "normstress.Set (Norm(model.stress))\n",
    "\n",
    "Draw (normstress, mesh, \"mises\")\n",
    "SetVisualization (min=0, max=20, deformation=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Navier Stokes model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "geo = netgen.geom2d.SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (2, 0.41), bcs = (\"wall\", \"outlet\", \"wall\", \"inlet\"))\n",
    "geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc=\"cyl\", maxh=0.02)\n",
    "mesh = Mesh( geo.GenerateMesh(maxh=0.07))\n",
    "Draw (mesh)\n",
    "mesh.Curve(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "timestep = 0.001\n",
    "navstokes = NavierStokes (mesh, nu=0.001, order=4, timestep = timestep,\n",
    "                              inflow=\"inlet\", outflow=\"outlet\", wall=\"wall|cyl\",\n",
    "                              uin=CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) ))\n",
    "                              \n",
    "\n",
    "navstokes.SolveInitial()\n",
    "\n",
    "Draw (navstokes.velocity, mesh, \"velocity\", reset=True)\n",
    "Draw (navstokes.pressure, mesh, \"pressure\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tend = 2\n",
    "t = 0\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        navstokes.DoTimeStep()\n",
    "        t = t+timestep\n",
    "        Redraw(blocking=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Transport equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh( netgen.geom2d.unit_square.GenerateMesh(maxh=0.1))\n",
    "\n",
    "timestep = 3e-3\n",
    "transport = TransportEquation (mesh, order=6, wind = (y-0.5, -x+0.5), timestep=timestep)\n",
    "\n",
    "transport.SetInitial( exp (-100* ((x-0.8)**2 + (y-0.5)**2) ) )\n",
    "\n",
    "Draw (transport.concentration, mesh, \"c\", reset=True)\n",
    "SetVisualization (min=0, max=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tend = 10\n",
    "t = 0\n",
    "with TaskManager(pajetrace=100*1000*1000):\n",
    "    while t < tend:\n",
    "        transport.DoTimeStep()\n",
    "        t = t+timestep\n",
    "        Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Time-dependent Convection-Diffusion Equation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh( netgen.geom2d.unit_square.GenerateMesh(maxh=0.1))\n",
    "\n",
    "timestep = 10e-3\n",
    "convdiff = ConvectionDiffusionEquation (mesh, order=3, lam=1e-3, wind = (y-0.5, -x+0.5), dirichlet=\".*\", timestep=timestep)\n",
    "\n",
    "convdiff.SetInitial( exp (-100* ((x-0.8)**2 + (y-0.5)**2) ) )\n",
    "Draw (convdiff.concentration, mesh, \"c\", reset=True)\n",
    "SetVisualization (min=0, max=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "tend = 10\n",
    "t = 0\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        convdiff.DoTimeStep()\n",
    "        t = t+timestep\n",
    "        Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Natural Convection:\n",
    "\n",
    "Coupling of Navier-Stokes and heat transport\n",
    "\n",
    "* Change in temperature leads to gravity forces\n",
    "* Temperature is convected by fluid velocity\n",
    "\n",
    "Rayleigh-Benard benchmark example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = netgen.geom2d.SplineGeometry()\n",
    "geo.AddRectangle( (0,0), (0.06, 0.01), bcs=['b','r','t','l'])\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.002))\n",
    "Draw (mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "timestep = 0.5\n",
    "navstokes = NavierStokes (mesh, nu=1.04177e-6, order=3, timestep = timestep,\n",
    "                              inflow=\"\", outflow=\"\", wall=\"l|r|b|t\", uin=(0,0) )\n",
    "\n",
    "Tinitial = 293.5-50*y+y*(0.01-y)*1e3*sin(20/0.06*x*math.pi)\n",
    "\n",
    "convdiff = ConvectionDiffusionEquation (mesh, order=3, lam=1.38e-7, \\\n",
    "                wind = navstokes.velocity, dirichlet=\"b|t\", udir=Tinitial, timestep=timestep)\n",
    "convdiff.SetInitial(Tinitial)\n",
    "\n",
    "T0 = 293\n",
    "beta = 2.07e-4\n",
    "navstokes.AddForce ( (1-beta*(convdiff.concentration-T0))*(0, -9.81))\n",
    "\n",
    "navstokes.SolveInitial()\n",
    "\n",
    "Draw (navstokes.pressure, mesh, \"pressure\", reset=True)\n",
    "Draw (navstokes.velocity, mesh, \"velocity\")\n",
    "visoptions.scalfunction='velocity:0'\n",
    "Draw (convdiff.concentration, mesh, \"temp\")\n",
    "input (\"key\")\n",
    "t, tend = 0, 1000\n",
    "with TaskManager(pajetrace=100*1000*1000):\n",
    "    while t < tend:\n",
    "        print (t)\n",
    "        navstokes.DoTimeStep()\n",
    "        convdiff.DoTimeStep()\n",
    "        t = t+timestep\n",
    "        Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Have a look into some of the templates:\n",
    "    \n",
    "   * Elasticity \n",
    "   * ConvectionDiffusion \n",
    "   * NavierStokesSIMPLE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Ongoing work:\n",
    "- unified interfaces to user and solvers\n",
    "- GUI integration  (--> Christopher Lackner)\n",
    "- everything with ALE (--> Michael Neunteufel)\n",
    "- same solvers in parallel  (--> Lukas Kogler)\n",
    "- higher order time-stepping / space time \n",
    "- shape optimization ( --> Peter Gangl + Kevin Sturm)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# How to contribute ? \n",
    "\n",
    "- Try NGSolve Model Templates\n",
    " - Suggest missing features \n",
    " - Suggest how to implement missing features\n",
    "- field experts needed:\n",
    " - Cooperate on research on numerical methods like\n",
    "   solvers, hpc, time-stepping, space-time, optimization, model order reduction\n",
    "  \n",
    " - Cooperate on research in applications\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# the end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}