{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.2 Constructive Solid Geometry (CSG)\n",
    "These geometries are bases on primitives (e.g. sphere, cylinder, plane) which are used to build solids by performing boolean operations. Netgen offers the following primitives\n",
    "\n",
    "| primitive  |  csg syntax |   meaning   |\n",
    "|:-----------|:------------|:------------|\n",
    "| half-space | Plane(Pnt a,Vec n)  |     point p in plane, normal vector    \n",
    "| sphere     | Sphere(Pnt c,float r)|    sphere with center c and radius r \n",
    "| cylinder   | Cylinder(Pnt a, Pnt b, float r) |  points a and b define the axes of a infinite cylinder of radius r \n",
    "| brick      | OrthoBrick (Pnt a, Pnt b) | axes parallel brick with minimal coordinates a and maximal coordinates b\n",
    "\n",
    "and the boolean operators\n",
    "\n",
    "| operator  |  set operation |\n",
    "|:----------|:---------------|\n",
    "| $*$ \t| intersection\n",
    "| $+$ \t| union\n",
    "| $-$ \t| intersection with complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import netgen.gui\n",
    "from ngsolve import Draw, Redraw # just for visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these primitives and operations, we can easily construct a cube. First we import the `netgen.csg` module, create 6 plane and intersect them to get the solid `cube`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import *\n",
    "\n",
    "left  = Plane (Pnt(0,0,0), Vec(-1,0,0) )\n",
    "right = Plane (Pnt(1,1,1), Vec( 1,0,0) )\n",
    "front = Plane (Pnt(0,0,0), Vec(0,-1,0) )\n",
    "back  = Plane (Pnt(1,1,1), Vec(0, 1,0) )\n",
    "bot   = Plane (Pnt(0,0,0), Vec(0,0,-1) )\n",
    "top   = Plane (Pnt(1,1,1), Vec(0,0, 1) )\n",
    "\n",
    "cube = left * right * front * back * bot * top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create a `CSGeometry` object and add the solid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSGeometry()\n",
    "geo.Add (cube)\n",
    "\n",
    "mesh = geo.GenerateMesh(maxh=0.25)\n",
    "Redraw()\n",
    "# mesh.Save(\"cube.vol\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import *\n",
    "\n",
    "cube = OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) )\n",
    "hole = Cylinder ( Pnt(0.5, 0.5, 0), Pnt(0.5, 0.5, 1), 0.2)\n",
    "\n",
    "geo = CSGeometry()\n",
    "geo.Add (cube-hole.maxh(0.05))\n",
    "mesh = geo.GenerateMesh(maxh=0.1)\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting properties of solids\n",
    "A solid has members which we can set to define the desired properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere = Sphere(Pnt(0,0,0),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can set a boundary name and a maximal mesh size on the surface of this sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere.bc(\"sphere\").maxh(0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define a material for the volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere.mat(\"iron\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case we want to visualize the geometry we can define the color (using rgb values) and transparency of the solid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere.col([1,0,0])#.transp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSGeometry()\n",
    "geo.Add(sphere)\n",
    "geo.Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = geo.GenerateMesh()\n",
    "print(type(ngmesh))\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To improve the approximation of curved geometries it is possible to use curved elements. This can be done within `NGSolve`. Thus we have to convert the `Netgen` mesh to a `NGSolve` mesh before curving it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.comp import Mesh\n",
    "mesh = Mesh(ngmesh)\n",
    "print(type(mesh))\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Curve(3)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the mesh size\n",
    "There are the following options to set the mesh size:\n",
    "\n",
    "* globally as argument `maxh` of `GenerateMesh`\n",
    "* to the surface of one solid (`maxh` property as above mentioned)\n",
    "* for the volume of a solid as optional argument when adding it to the geometry `Add(...,bc)`\n",
    "* restrict the mesh size for one point using `RestrictH`\n",
    "* use `CloseSurfaces` to generate anisotropic meshes\n",
    "\n",
    "### Global mesh size\n",
    "The global mesh size can be set with the named argument `maxh`. The following two versions are equivalent since all arguments of the of the `GenerateMesh` function are parsed to the `MeshingParameters` if no named argument `mp` is given."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit_cube.GenerateMesh(maxh=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.meshing import MeshingParameters\n",
    "mp = MeshingParameters(maxh=0.4)\n",
    "unit_cube.GenerateMesh(mp = mp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh size for one solid\n",
    "To set the mesh size for one domain of the mesh we have to add the desired `maxh` as argument when adding the solid to the geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSGeometry()\n",
    "\n",
    "brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))\n",
    "sphere = Sphere(Pnt(0,0,0),1)\n",
    "\n",
    "geo.Add(brick-sphere)\n",
    "geo.Add(sphere,maxh=0.1)\n",
    "ngmesh = geo.GenerateMesh(maxh=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh size on a surface\n",
    "If we want to refine just on a surface we define it as property of the solid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSGeometry()\n",
    "\n",
    "brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))\n",
    "sphere = Sphere(Pnt(0,0,0),1)\n",
    "\n",
    "geo.Add(brick-sphere)\n",
    "geo.Add(sphere.maxh(0.1))\n",
    "ngmesh = geo.GenerateMesh()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh size in points\n",
    "This can be done with the `MeshingParameters`. Using `RestrictH` we can define the mesh size in an arbitrary point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSGeometry()\n",
    "\n",
    "brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))\n",
    "sphere = Sphere(Pnt(0,0,0),1)\n",
    "\n",
    "mp = MeshingParameters(maxh=0.4)\n",
    "mp.RestrictH (x=0, y=0, z=1, h=0.025)\n",
    "        \n",
    "geo.Add(brick-sphere)\n",
    "geo.Add(sphere)\n",
    "ngmesh = geo.GenerateMesh(mp = mp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Anisotropic meshes\n",
    "If the geometry contains thin layers we can use `CloseSurfaces` to avoid elements with small angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import *\n",
    "\n",
    "geo = CSGeometry()\n",
    "\n",
    "box = OrthoBrick(Pnt(0,0,0),Pnt(1,1,1))\n",
    "top = Plane(Pnt(0,0,0.52),Vec(0,0,1))\n",
    "bot = Plane(Pnt(0,0,0.48),Vec(0,0,-1))\n",
    "plate = box * top * bot\n",
    "\n",
    "geo.Add((box-top).mat(\"air\"))\n",
    "geo.Add(plate.mat(\"plate\"))\n",
    "geo.Add((box-bot).mat(\"air\"))\n",
    "\n",
    "slices = [2**(-i) for i in reversed(range(1,6))]\n",
    "# define the close surfaces\n",
    "geo.CloseSurfaces(bot,top)#,slices)\n",
    "nmesh = geo.GenerateMesh(maxh=0.3)\n",
    "# refine the mesh between the close surfaces\n",
    "# ZRefinement(nmesh,geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting boundary conditions\n",
    "### Boundary condition on the surface of a solid\n",
    "Setting a boundary condition on the whole surface of solid can be achieved by adding it as property to the solid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2)).bc('outer')\n",
    "sphere = Sphere(Pnt(0,0,0),1).bc('sphere')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modify boundary between two solids\n",
    "This can be done by adding the named argument `bcmod` when adding the solid to the geometry. Here we change the boundary condition on the surface between the `halfsphere` and the already added `box`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "halfsphere = sphere * Plane(Pnt(0,0,0),Vec(1,0,0)).bc('plane')\n",
    "box = brick-sphere\n",
    "geo = CSGeometry()\n",
    "geo.Add(box.col([1,0,0]).transp())\n",
    "geo.Add(halfsphere.col([0,0,1]),bcmod=[(box,\"halfsphere\")])\n",
    "geo.Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = geo.GenerateMesh()\n",
    "mesh = Mesh(ngmesh)\n",
    "mesh.GetBoundaries()"
   ]
  },
  {
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}