2D_Periodic_Mesh_with_GMSH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "38235d78",
   "metadata": {},
   "source": [
    "# Periodic Mesh in GMSH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8d65fc69",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gmsh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "163dcf60",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ca26ec2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.add(\"Periodic x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "53eb03ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "minCoords = (0., 0.)\n",
    "maxCoords = (1.0, 1.0)\n",
    "cellSize = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f7ec766c",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundaries = {\n",
    "   \"Bottom\": 1,\n",
    "   \"Top\": 2,\n",
    "   \"Right\": 3,\n",
    "   \"Left\": 4,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "39ff9923",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, ymin = minCoords\n",
    "xmax, ymax = maxCoords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ad5595f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = gmsh.model.geo.add_point(xmin,ymin,0., meshSize=cellSize)\n",
    "p2 = gmsh.model.geo.add_point(xmax,ymin,0., meshSize=cellSize)\n",
    "p3 = gmsh.model.geo.add_point(xmin,ymax,0., meshSize=cellSize)\n",
    "p4 = gmsh.model.geo.add_point(xmax,ymax,0., meshSize=cellSize)\n",
    "\n",
    "l1 = gmsh.model.geo.add_line(p1, p2, tag=boundaries[\"Bottom\"])\n",
    "l2 = gmsh.model.geo.add_line(p2, p4, tag=boundaries[\"Right\"])\n",
    "l3 = gmsh.model.geo.add_line(p4, p3, tag=boundaries[\"Top\"])\n",
    "l4 = gmsh.model.geo.add_line(p3, p1, tag=boundaries[\"Left\"])\n",
    "\n",
    "cl = gmsh.model.geo.add_curve_loop((l1, l2, l3, l4))\n",
    "surface = gmsh.model.geo.add_plane_surface([cl])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "55655c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.geo.synchronize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aec190f3",
   "metadata": {},
   "source": [
    "The periodicity transform is provided as a 4x4 affine transformation matrix, given by row.\n",
    "During mesh generation, the mesh on the right side will be created by copying the mesh from the left side.\n",
    "\n",
    "To impose that the right side should match the left side, the following periodicity constraint is set:\n",
    "\n",
    "$M = \\begin{pmatrix}\n",
    "  1 & 0 & 0 & 1\\\\ \n",
    "  0 & 1 & 0 & 0\\\\\n",
    "  0 & 0 & 1 & 0\\\\\n",
    "  0 & 0 & 0 & 1\n",
    "\\end{pmatrix}$\n",
    "\n",
    "The left side is copied by translating it along the x axis.\n",
    "\n",
    "Such a 4 by 4 matrix M corresponds to a afﬁne transformation T() that transforms point (or vector) x to point (or vector) y. The upper-left 3 × 3 sub-matrix of the matrix shown above represents a rotation transform, but may also include scales and shears. The last column of the matrix represents a translation. When used as a coordinate system, the upper-left 3 x 3 sub-matrix represents an orientation in space while the last column vector represents a position in space. The transformation T() of point x to point y is obtained by performing the matrix-vector multiplication:\n",
    "\n",
    "$ y = Mx $ \n",
    "\n",
    "see for more information:\n",
    "\n",
    "https://www.brainvoyager.com/bv/doc/UsersGuide/CoordsAndTransforms/SpatialTransformationMatrices.html#:~:text=The%204%20by%204%20transformation,in%20the%20first%20three%20columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8ce7e1d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "translation = [1, 0, 0, 1, \n",
    "               0, 1, 0, 0,\n",
    "               0, 0, 1, 0, \n",
    "               0, 0, 0, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8370444a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.mesh.setPeriodic(1, [boundaries[\"Right\"]], [boundaries[\"Left\"]], translation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "99de8256",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add Physical groups\n",
    "for name, tag in boundaries.items():\n",
    "    gmsh.model.add_physical_group(1, [tag] , tag)\n",
    "    gmsh.model.set_physical_name(1, tag, name)\n",
    "\n",
    "gmsh.model.addPhysicalGroup(2, [surface], surface)\n",
    "gmsh.model.setPhysicalName(2, surface, \"Elements\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c97212b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4], dtype=int32)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmsh.model.mesh.get_periodic(1, [boundaries[\"Right\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e0d87851",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4], dtype=int32)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmsh.model.mesh.get_periodic(1, [boundaries[\"Left\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3bbdf1bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#gmsh.fltk.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f9354762",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : Reconstructing periodicity for curve connection 3 - 4\n",
      "Info    : Meshing 1D...\n",
      "Info    : [  0%] Meshing curve 1 (Line)\n",
      "Info    : [ 30%] Meshing curve 2 (Line)\n",
      "Info    : [ 80%] Meshing curve 4 (Line)\n",
      "Info    : [100%] Meshing curve 3 (Line) as a copy of curve 4\n",
      "Info    : Done meshing 1D (Wall 0.00341694s, CPU 0.005503s)\n",
      "Info    : Meshing 2D...\n",
      "Info    : Meshing surface 1 (Plane, Frontal-Delaunay)\n",
      "Info    : Done meshing 2D (Wall 0.022048s, CPU 0.019632s)\n",
      "Info    : Reconstructing periodicity for curve connection 3 - 4\n",
      "Info    : 142 nodes 286 elements\n",
      "Info    : Writing 'periodicx.msh'...\n",
      "Info    : Done writing 'periodicx.msh'\n"
     ]
    }
   ],
   "source": [
    "gmsh.model.mesh.generate(2) \n",
    "gmsh.write(\"periodicx.msh\")\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75bee851",
   "metadata": {},
   "source": [
    "## Importing in PETSc\n",
    "\n",
    "Could be useful:\n",
    "https://github.com/firedrakeproject/firedrake/issues/1852"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "26fd2c5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petsc4py import PETSc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "3df14bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plex = PETSc.DMPlex().createFromFile(\"periodicx.msh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "94df03ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DM Object: unnamed 1 MPI process\n",
      "  type: plex\n",
      "unnamed in 2 dimensions:\n",
      "  Number of 0-cells per rank: 131\n",
      "  Number of 1-cells per rank: 373\n",
      "  Number of 2-cells per rank: 242\n",
      "Periodic mesh coordinates localized\n",
      "Labels:\n",
      "  celltype: 3 strata with value/size (0 (131), 3 (242), 1 (373))\n",
      "  depth: 3 strata with value/size (0 (131), 1 (373), 2 (242))\n",
      "  Cell Sets: 1 strata with value/size (1 (242))\n",
      "  Face Sets: 4 strata with value/size (1 (10), 2 (10), 3 (10), 4 (10))\n"
     ]
    }
   ],
   "source": [
    "plex.view()"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}