{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.3 Working with meshes\n",
    "## One-dimensional meshes\n",
    "Meshes in one-dimension can be constructed using the `netgen.meshing` module.\n",
    "We just have to add segments (`Element1D`) and the boundaries (`Element0D`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import netgen.gui\n",
    "from netgen.meshing import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create a new `Mesh` and set the spatial dimension to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Mesh(dim=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define and add the `MeshPoint`'s to the mesh. The function `m.Add` returns `PointId`'s which we store in an array to be able to construct the segments in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "pnums = []\n",
    "for i in range(0, N+1):\n",
    "    pnums.append (m.Add (MeshPoint (Pnt(i/N, 0, 0))))\n",
    "\n",
    "type(pnums[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can loop over the array with the `PointId`'s and add one-dimensional elements to the mesh. Further we can set the material for our domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = m.AddRegion(\"material\", dim=1)\n",
    "for i in range(0,N):\n",
    "    m.Add (Element1D ([pnums[i],pnums[i+1]], index=idx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we have to add the boundary elements and set boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_left = m.AddRegion(\"left\", dim=0)\n",
    "idx_right = m.AddRegion(\"right\", dim=0)\n",
    "\n",
    "m.Add (Element0D (pnums[0], index=idx_left))\n",
    "m.Add (Element0D (pnums[N], index=idx_right))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To be able to visualize one-dimensional meshes and solution activate `Show edges` in the menu `View > Viewing options > Mesh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "mesh = ngsolve.Mesh(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-dimensional meshes\n",
    "As example we mesh a unit square [0,1]x[0,1] using quadrilaterals.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create an empty mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = Mesh(dim=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and add all the `MeshPoint`'s we will need for the final mesh. Similar to the one-dimensional mesh we store the `PointId`'s in the `pnums` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=5\n",
    "pnums = []\n",
    "for i in range(N + 1):\n",
    "    for j in range(N + 1):\n",
    "        pnums.append(ngmesh.Add(MeshPoint(Pnt(i / N, j / N, 0))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a region, and add the quadrilaterals to the mesh. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_dom = ngmesh.AddRegion(\"mat\", dim=2)\n",
    "for j in range(N):\n",
    "    for i in range(N):\n",
    "        ngmesh.Add(Element2D(idx_dom, [pnums[i + j * (N + 1)],\n",
    "                               pnums[i + (j + 1) * (N + 1)],\n",
    "                               pnums[i + 1 + (j + 1) * (N + 1)],\n",
    "                               pnums[i + 1 + j * (N + 1)]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we have to add boundary elements and set boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# horizontal boundaries\n",
    "for i in range(N):\n",
    "   ngmesh.Add(Element1D([pnums[N + i * (N + 1)],\n",
    "                       pnums[N + (i + 1) * (N + 1)]], index=1))\n",
    "   ngmesh.Add(Element1D([pnums[0 + i * (N + 1)], pnums[0 + (i + 1) * (N + 1)]], index=1))\n",
    "\n",
    "# vertical boundaries\n",
    "for i in range(N):\n",
    "   ngmesh.Add(Element1D([pnums[i], pnums[i + 1]], index=2))\n",
    "   ngmesh.Add(Element1D([pnums[i + N * (N + 1)], pnums[i + 1 + N * (N + 1)]], index=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.webgui import Draw\n",
    "mesh = ngsolve.Mesh(ngmesh)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge three-dimensional meshes\n",
    "\n",
    "In the following example we will merge two surface meshes and generate a unified volume mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.meshing import *\n",
    "from netgen.csg import *\n",
    "\n",
    "from ngsolve import ngsglobals\n",
    "ngsglobals.msg_level = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As starting point we create two geometries and mesh them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate brick and mesh it\n",
    "geo1 = CSGeometry()\n",
    "geo1.Add (OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) ))\n",
    "m1 = geo1.GenerateMesh (maxh=0.1)\n",
    "# m1.Refine()\n",
    "\n",
    "# generate sphere and mesh it\n",
    "geo2 = CSGeometry()\n",
    "geo2.Add (Sphere (Pnt(0.5,0.5,0.5), 0.1))\n",
    "m2 = geo2.GenerateMesh (maxh=0.05)\n",
    "m2.Refine()\n",
    "# m2.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we start the merging process. Therefore we create an empty mesh and add a `FaceDescriptor` for each of the surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an empty mesh\n",
    "ngmesh = Mesh()\n",
    "\n",
    "# a face-descriptor stores properties associated with a set of surface elements\n",
    "# bc .. boundary condition marker,\n",
    "# domin/domout .. domain-number in front/back of surface elements (0 = void),\n",
    "# surfnr .. number of the surface described by the face-descriptor\n",
    "\n",
    "fd_outside = ngmesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1))\n",
    "fd_inside = ngmesh.Add (FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the surface elements stay the same in the merged mesh, we copy the points on the surface and the surface elements to the new mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy all boundary points from first mesh to new mesh.\n",
    "# pmap1 maps point-numbers from old to new mesh\n",
    "pmap1 = { }\n",
    "for e in m1.Elements2D():\n",
    "    for v in e.vertices:\n",
    "        if (v not in pmap1):\n",
    "            pmap1[v] = ngmesh.Add (m1[v])\n",
    "\n",
    "\n",
    "# copy surface elements from first mesh to new mesh\n",
    "# we have to map point-numbers:\n",
    "\n",
    "for e in m1.Elements2D():\n",
    "    ngmesh.Add (Element2D (fd_outside, [pmap1[v] for v in e.vertices]))\n",
    "\n",
    "# same for the second mesh:\n",
    "pmap2 = { }\n",
    "for e in m2.Elements2D():\n",
    "    for v in e.vertices:\n",
    "        if (v not in pmap2):\n",
    "            pmap2[v] = ngmesh.Add (m2[v])\n",
    "\n",
    "for e in m2.Elements2D():\n",
    "    ngmesh.Add (Element2D (fd_inside, [pmap2[v] for v in e.vertices]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we have to generate the new volume mesh. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh.GenerateVolumeMesh()\n",
    "import ngsolve\n",
    "mesh = ngsolve.Mesh(ngmesh)\n",
    "Draw(mesh);"
   ]
  },
  {
   "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
}