NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "(basics:modelling-domains-with-occ)=\n",
    "# Open Cascade Technology Geometry (OCC Geometry)\n",
    "\n",
    "Netgen provides a Python wrapper around the Open Cascade Technology (OCCT) geometry kernel. It allows to model complex geometric objects. It also allows to import models via STEP format, explore, and modify the geometry. \n",
    "\n",
    "In case you are building Netgen/NGsolve from scratch you need the cmake-flag `-DUSE_OCC=ON`.\n",
    "\n",
    "In the wrapper we aim in bringing most of the structure of OCCT to Python. If you are familiar with the C++ interface of Open Cascade you will recognize many classes. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3febe816",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a253516",
   "metadata": {},
   "source": [
    "## Construction of 2D objects (Workplanes)\n",
    "\n",
    "Workplanes are the basic building blocks of 2D geometry, they can be understood as a plane in 3D space.\n",
    "\n",
    "There are several methods that are already implemented to define standard shapes (`Rectangle`, `Circle`, `Ellipse`...).\n",
    "\n",
    "After \"sketching\" on the workplane you need to extact the object of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8ee20622",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<netgen.libngpy._NgOCC.WorkPlane at 0x183ace664f0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wp = WorkPlane()\n",
    "wp.Rectangle(2, 5) # you still cannot draw a rectangle with the same syntax as in the GUI\n",
    "\n",
    "# Draw(topo_rect) # try ... it fails"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a5cd609",
   "metadata": {},
   "source": [
    "Before drawing the object you need to extract the face from the workplane. This is done by the `Face` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bc4cc832",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "55ac3a67b92f4ed7b555ee800bd41f14",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rect = wp.Face()\n",
    "Draw(rect);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16dab61",
   "metadata": {},
   "source": [
    "### Orientation of the workplane\n",
    "\n",
    "The concept of Workplane is quite powerful, it allows to draw on a 2D plane and automatically map it to the 3D space.\n",
    "\n",
    "The `WorkPlane` is initialized with:\n",
    "\n",
    "```python\n",
    "class WorkPlane(\n",
    "    axes: Axes = ...,\n",
    "    pos: gp_Ax2d = ...\n",
    ")\n",
    "```\n",
    "\n",
    "where `Axes` is a set a point and two vectors that define the orientation of the workplane. \n",
    "\n",
    "```python\n",
    "Axes(\n",
    "    p: gp_Pnt = (0, 0, 0), # origin \n",
    "    n: gp_Dir = (0, 0, 1), # normal vector\n",
    "    h: gp_Dir = (1, 0, 0)  # horizontal vector\n",
    "    ) -> None\n",
    "```\n",
    "\n",
    "The `gp_Ax2d` is a set of a point in 2D and the main direction.\n",
    "\n",
    "```python\n",
    "gp_Ax2d(\n",
    "    p: gp_Pnt2d = (0, 0), # origin\n",
    "    v: gp_Dir2d = (1, 0)  # direction\n",
    "    ) -> None\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c3b6b79a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3d695b5adeea4635811ed195f4af2a31",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axes = Axes( (1,1,0) ,(0,1,0) , (-1,1,0) ) \n",
    "\n",
    "# axes = Axes((1, 1, 0), Y, Y-X) # also works\n",
    "\n",
    "pos = gp_Ax2d()\n",
    "\n",
    "wp2 = WorkPlane( axes = axes, pos = pos )\n",
    "rect2 = wp2.Rectangle(2, 5).Face()\n",
    "Draw(rect2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f520502b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5fc67deb767142998079c577187f1760",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "union = rect + rect2\n",
    "Draw(union);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06d4c3b0",
   "metadata": {},
   "source": [
    "\n",
    "Let's create the a workplane, the default position is the center of the plane and the default direction is the x-axis.\n",
    "\n",
    "To create a figure think of it like giving directions to a driver. For example: \n",
    "\n",
    "    1. Go straight for 1 meter\n",
    "    2. immediately turn (left) 60 degrees\n",
    "\n",
    "then repeat the previous steps for 6 times in total. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1fcc5e74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1fc7e1cc8ea4fbca86fc150aebdfcd7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wp = WorkPlane()\n",
    "for i in range(6):\n",
    "    wp.Line(1).Rotate(60)\n",
    "Draw(wp.Face());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79435494",
   "metadata": {},
   "source": [
    "Suppose to change slightly the above instructions to:\n",
    "\n",
    "    1. Go straight for 1 meter\n",
    "    2. Turn (right) 60 degrees while still driving so that you make an arc of 0.5 meters\n",
    "\n",
    "then repeat the previous steps for 6 times in total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e72fdb16",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f32f660118d646daac4edc0f0e528a16",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(6):\n",
    "    wp.Line(1).Arc(0.5, 60)\n",
    "face = wp.Face()\n",
    "Draw(face);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b47f90c",
   "metadata": {},
   "source": [
    "The drawing needs to be positively oriented, otherwise the face will not be created.\n",
    "\n",
    "If you have drawn the object in the wrong direction you can use the `Reverse` method to change the orientation.\n",
    "\n",
    "A negative oriented Wire is the same of subtracting to the figure a positive oriented Wire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6fe9a313",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "82cde27ef9b344f095930459815df207",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wp = WorkPlane().RectangleC(2, 1) \\\n",
    "    .Circle(0.5, 0, 0.2).Reverse() \\\n",
    "    .Circle(-0.5, 0, 0.2).Reverse()\n",
    "Draw(wp.Face());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f541961",
   "metadata": {},
   "source": [
    "There are two main ways to obtain a 3D object from a 2D object: \n",
    "\n",
    "1. Extrusion in a certain direction\n",
    "2. Revolution around an axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6b3e417b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d52235289db4c2eb8f090e597d52709",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "face = WorkPlane().MoveTo(10, 0).Line(18).Arc(2, 90).Line(6).Arc(2, 90) \\\n",
    "    .Line(18).Rotate(90).Line(2).Rotate(90).Line(10) \\\n",
    "    .Arc(3, -180).Line(10, \"nice line\").Close().Face()\n",
    "\n",
    "Draw(face);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1ecf8a49",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "de415c58d660405ca0d6d08ae81faeee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "ext = face.Extrude(-40).Move((10, 0, 0))\n",
    "\n",
    "rev = face.Revolve(Axis((0, 0, 0), Y+Z+X), 280).Move((-10, 0, 0))\n",
    "\n",
    "Draw(rev+ext );"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46853e07",
   "metadata": {},
   "source": [
    "### Modify 3D objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ad607800",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "697fb34083464f7691d88ca00465591c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fillet = ext.MakeFillet(ext.edges, 0.3)\n",
    "\n",
    "Draw(fillet);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d22a824d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "869b5113f817406d9d29ac03e99b9afd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "chamfer = ext.MakeChamfer(ext.edges, 0.3)\n",
    "\n",
    "Draw(chamfer);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad552dc0",
   "metadata": {},
   "source": [
    "There are several other ready-to-use 2D shapes. Look into the module to find more."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "### Construction of 3D objects\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "We define a box aligned with the Cartesian coordinates given by two points with minimal and maximal x/y/z coordinates.  A cylinder is given by a point on the axis, a direction vector, the radius $r$, and the height $h$. The symbols `X`, `Y`, and `Z` are predefined basis vectors for the Cartesian coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "box = Box(Pnt(0,0,0), Pnt(1,1,1))\n",
    "cyl = Cylinder(Pnt(1,0.5,0.5), X, r=0.3, h=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "The Boolean operations `fuse`, `common`, and `cut` provided by OCCT are made available by the operators `+`, `*` and `-`. More can be found in the [OCCT-documentation](https://dev.opencascade.org/doc/refman/html/class_b_rep_algo_a_p_i___boolean_operation.html). Note that the fuse generates one new solid, there is no interface face where the cylinder is touching the box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a2355d099624f80887c1a52122d8553",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "# name the faces \n",
    "cyl.faces.Min(X).name = \"cyl_minx\"\n",
    "box.faces.Min(X).name = \"box_minx\"\n",
    "box.faces.Max(X).name = \"box_maxx\"\n",
    "\n",
    "# create the union\n",
    "fused = box+cyl\n",
    "\n",
    "# impose beforehand the mesh size\n",
    "fused.faces[\"box_.*\"].maxh = 0.05\n",
    "\n",
    "# che check the names\n",
    "Draw(fused);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "The generated objects are Py-wrapped OCCT objects derived from TopoDS_Shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8bd478b4cf7d4dfd9ed6cf068066246e",
       "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"
    }
   ],
   "source": [
    "from ngsolve import Mesh\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "geo = OCCGeometry(fused)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.2, grading=0.9))\n",
    "mesh.Curve(4)\n",
    "Draw (mesh, clipping={\"z\":-1});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "Instead of fusing, we can glue together shapes. Then, the resulting composite solid contains the interface face between the solids. This is important when dealing with separate material regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c549788d762c490ba17db1689ce42bba",
       "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"
    }
   ],
   "source": [
    "geo = Glue( [ cyl, box]) \n",
    "\n",
    "mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.2, grading=0.9))\n",
    "mesh.Curve(4)\n",
    "\n",
    "Draw (mesh, clipping={\"z\": -1});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab49259",
   "metadata": {},
   "source": [
    "**be aware of the order**\\\n",
    "change the order and look at the names of the common face"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "A third option is to form a compound of shapes, then the component shapes are meshed independently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f78975ceb7c4b9daf7460af8847d624",
       "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"
    }
   ],
   "source": [
    "geo = Compound( [box, cyl])\n",
    "\n",
    "mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.2)).Curve(3)\n",
    "Draw (mesh, clipping=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Transformation of shapes\n",
    "\n",
    "We can translate, rotate and mirror shapes. The transformation preserves the original shape, and returns a transformed copy. The translation `Move` takes a vector as argument, the rotation `Rotation` needs an axis given by a point and a direction, and an angle. The sign of the angle reflects the right hand screw rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "16",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "af8b1bee56ed4f9daef884688544bb2d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solid = Box((0,0,0), (5,3,1)) + Sphere((0,0,0), 0.3)\n",
    "solid2 = solid.Move((5,0,2))\n",
    "solid3 = solid.Move((0,0,4)).Rotate( Axis((0,0,4), X), 45)\n",
    "Draw (solid + solid2 + solid3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## Selection of shapes\n",
    "\n",
    "Sometimes we want to set properties of shapes, boundary conditions or mesh-size. For that we can use shape selectors. \n",
    "\n",
    "* The `Max` and `Min` selectors finds the sub-shapes where the center of gravity has maximal or minimal coordinates in a given direction\n",
    "\n",
    "More of them will come soon ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f8c90bb981444854b900c74fedc49dc1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "box = Box((0,0,0), (1,1,1))\n",
    "box.faces.Max(Z).col = (1,0,0)\n",
    "box.faces.Min(Y).col = (0,0,1)\n",
    "Draw(box);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "Faces where center of gravity has $x$ coordinate less than 0.8:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e2919f7a8afe40b58c742ad0c22e8168",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'ngsolve_version': 'Netgen x.x', 'mesh_dim': …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Draw(Compound(box.faces[X<0.8]));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "### Last remark on the workflow OCC-Netgen-NGSolve\n",
    "\n",
    "The workflow on how to create a domain is as follows:\n",
    "1. Create the geometry using OCC \n",
    "1. Pass the geometry to Netgen with the meshing parameters\n",
    "2. Mesh the geometry using Netgen \n",
    "3. Pass the netgen mesh to NGSolve \n",
    "\n",
    "\n",
    "| OCC shape| Netgen shape| Netgen mesh | NGSolve mesh |\n",
    "| --- | ------ | ------- | --- |\n",
    "| TopoDS\\_Shape | OCCGeometry(...) | .GenerateMesh(...) |Mesh(...) |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7bdd731",
   "metadata": {},
   "source": [
    "## Metamaterials:\n",
    "\n",
    "using the OCC module it is possible to create complex geometries for metamaterials. Let's createa simle geometry with a certain chirality.\n",
    "\n",
    "in this case a crown with some additional arms that pop out of the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f32b7811",
   "metadata": {},
   "outputs": [],
   "source": [
    "def BaseShape(r_in, r_out, l):\n",
    "    c_in = Circle(Pnt(0, 0), r_in).Face()\n",
    "    c_out = Circle(Pnt(0, 0), r_out).Face()\n",
    "    h = r_out - r_in\n",
    "\n",
    "    r_top = MoveTo(0, r_in).Rectangle(l/2, h).Face()\n",
    "    t_right = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 90)\n",
    "    r_bot = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 180)\n",
    "    t_left = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 270)\n",
    "    shape = c_out - c_in + r_top + r_bot + t_right + t_left\n",
    "\n",
    "    return shape\n",
    "\n",
    "\n",
    "l = 6\n",
    "r_in = 1\n",
    "r_out = 2\n",
    "base = BaseShape(r_in, r_out, l)\n",
    "\n",
    "Draw(base);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f00b1eb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "figure = base\n",
    "for i in range(10):\n",
    "    for j in range(5):\n",
    "        if (i+j) % 2 == 0:\n",
    "            figure += base.Move((i*6, j*6, 0))\n",
    "        else:\n",
    "            figure += base.Mirror(Axis((0, 0, 0), (1, 0, 0))\n",
    "                                  ).Move((i*6, j*6, 0))\n",
    "\n",
    "\n",
    "Draw(figure)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f981caf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "figure = figure.Extrude(2)\n",
    "\n",
    "# add a left handle to the figure 1 x 5*l x 2\n",
    "lbar = Box((0, 0, 0), (1, 5*l, 2)).Move((-l/2-1, -l/2, 0))\n",
    "figure += lbar\n",
    "\n",
    "# add a right handle to the figure 1 x 5*l x 2\n",
    "rbar = Box((0, 0, 0), (1, 5*l, 2)).Move((9*l + l/2, -l/2, 0))\n",
    "figure += rbar\n",
    "\n",
    "figure.faces.Min(X).name = \"left\"\n",
    "figure.faces.Max(X).name = \"right\"\n",
    "\n",
    "\n",
    "# Draw(base)\n",
    "Draw(figure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3981935",
   "metadata": {},
   "source": [
    "**Example of elasticity problem with metamaterials**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5894fe1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.solvers import Newton\n",
    "\n",
    "\n",
    "E, nu = 210, 0.2\n",
    "mu = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "\n",
    "def C(u):\n",
    "    F = Id(3) + Grad(u)\n",
    "    return F.trans * F\n",
    "\n",
    "\n",
    "def NeoHooke(C):\n",
    "    return 0.5*mu*(Trace(C-Id(3)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)\n",
    "\n",
    "\n",
    "factor = Parameter(0)\n",
    "force = CoefficientFunction((factor, 0, 0))\n",
    "\n",
    "\n",
    "mesh = Mesh(OCCGeometry(figure).GenerateMesh(maxh=2)).Curve(4)\n",
    "\n",
    "fes = H1(mesh, order=2, dirichlet=\"left\", dim=mesh.dim)\n",
    "u = fes.TrialFunction()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += Variation(NeoHooke(C(u)).Compile()*dx)\n",
    "a += Variation((-InnerProduct(force, u)).Compile()*ds(\"right\"))\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:] = 0\n",
    "\n",
    "\n",
    "scene = Draw(gfu, mesh, deformation=True, scale=2)\n",
    "\n",
    "\n",
    "numsteps = 2\n",
    "for step in range(numsteps):\n",
    "    print(step)\n",
    "    factor.Set(0.25*(step+1)/numsteps)\n",
    "    with TaskManager():\n",
    "        Newton(a, gfu, printing=True, dampfactor=0.9)\n",
    "    scene.Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97c2536c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}