{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "18d98400-3dc3-4d11-9ed7-5197e0ba4c63",
   "metadata": {},
   "source": [
    "# [Gmsh 示例]{.lang-zh} [Gmsh Examples]{.lang-en}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} Overview\n",
    ":class: tip lang-en\n",
    "Generating meshes with the Gmsh Python API: a Möbius strip, extrusion for PML layers, and notes on Gmsh node ordering for high-order elements. Requires the `gmsh` package (this notebook is not executed during the book build).\n",
    "```"
   ],
   "id": "13f7d650"
  },
  {
   "cell_type": "markdown",
   "id": "17d899f7-3842-43db-af08-8ea000ca4e37",
   "metadata": {},
   "source": [
    "## Mobius strip"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2e259b5-2c9e-460c-a2d6-cbd0e806edce",
   "metadata": {},
   "source": [
    "### Plot by Matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a3747af-5b39-40b3-acf3-66650ef0de93",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gmsh\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# for interactive figures\n",
    "%matplotlib widget\n",
    "\n",
    "R = 10\n",
    "W = 5\n",
    "\n",
    "_u = np.linspace(0, 2*np.pi, 201, endpoint=True)\n",
    "_v = np.linspace(-1/2, 1/2, 10)\n",
    "u, v = np.meshgrid(_u, _v)\n",
    "\n",
    "r = R - W*v*np.sin(u/2)\n",
    "x = r*np.cos(u)\n",
    "y = r*np.sin(u)\n",
    "z = W*v*np.cos(u/2)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.axes(projection='3d')\n",
    "ax.plot_surface(x, y, z)\n",
    "zlim = ax.set_zlim([-2*W, 2*W])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad980a64-cf5a-48a4-b9b7-c42ff0b13c4d",
   "metadata": {},
   "source": [
    "### Python script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68f0ef52-90de-4af7-a023-eb5bfb7d0504",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import gmsh\n",
    "import numpy as np\n",
    "\n",
    "gmsh.initialize()\n",
    "fac = gmsh.model.geo\n",
    "R = 5\n",
    "W = 2\n",
    "M = 50\n",
    "_u = np.linspace(0, 2*np.pi, M+1, endpoint=True)\n",
    "_v = np.linspace(-1/2, 1/2, 2)\n",
    "u, v = np.meshgrid(_u, _v)\n",
    "\n",
    "r = R - W*v*np.sin(u/2)\n",
    "x = r*np.cos(u)\n",
    "y = r*np.sin(u)\n",
    "z = W*v*np.cos(u/2)\n",
    "\n",
    "points = [[], []]\n",
    "\n",
    "line_current = None\n",
    "ss = []\n",
    "for i in range(M):\n",
    "    points[0].append(fac.add_point(x[0, i], y[0, i], z[0, i]))\n",
    "    points[1].append(fac.add_point(x[1, i], y[1, i], z[1, i]))\n",
    "    line_pre = line_current\n",
    "    line_current = fac.add_line(points[0][-1], points[1][-1])\n",
    "    \n",
    "    if i == 0:\n",
    "        line0 = line_current\n",
    "    else:\n",
    "        j = i - 1\n",
    "        line3 = fac.add_line(points[0][i-1], points[0][i])\n",
    "        line4 = fac.add_line(points[1][i-1], points[1][i])\n",
    "\n",
    "\n",
    "        cl = fac.add_curve_loop([line_pre, line4, -line_current, -line3])\n",
    "        ss.append(fac.add_plane_surface([cl]))\n",
    "\n",
    "line3 = fac.add_line(points[0][-1], points[1][0])\n",
    "line4 = fac.add_line(points[1][-1], points[0][0])\n",
    "cl = fac.add_curve_loop([line_current, line4, line0, -line3])\n",
    "ss.append(fac.add_plane_surface([cl]))\n",
    "\n",
    "fac.add_surface_loop(ss)\n",
    "fac.synchronize()\n",
    "\n",
    "h = 0.4\n",
    "gmsh.option.setNumber(\"Mesh.MeshSizeMin\", h)\n",
    "gmsh.option.setNumber(\"Mesh.MeshSizeMax\", h)\n",
    "\n",
    "__old_verbosity = gmsh.option.getNumber(\"General.Verbosity\")\n",
    "gmsh.option.setNumber(\"General.Verbosity\", 1)\n",
    "# TODO: catch the exception of gmsh\n",
    "gmsh.model.mesh.generate()\n",
    "gmsh.option.setNumber(\"General.Verbosity\", __old_verbosity)\n",
    "\n",
    "# gmsh.write('gmsh/mobius.msh')\n",
    "\n",
    "# gui or not\n",
    "if '-popup' in sys.argv:\n",
    "    gmsh.fltk.run()\n",
    "\n",
    "\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08d7b607-4089-479a-90e8-9dba979a0caf",
   "metadata": {},
   "source": [
    "### Solve Poisson problem on Mobius strip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15a12f6a-7acd-4ebe-9c97-14b2fb216239",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from firedrake import *\n",
    "\n",
    "msh = Mesh('gmsh/mobius.msh', dim=3)\n",
    "triplot(msh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a09a5b1-4150-4d2c-8377-ef1fd0d73c99",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = SpatialCoordinate(msh)\n",
    "f = x**2 + y**2 + z**2\n",
    "\n",
    "V = FunctionSpace(msh, 'CG', 1)\n",
    "u, v = TrialFunction(V), TestFunction(V)\n",
    "a = inner(grad(u), grad(v))*dx - inner(f, v)*dx\n",
    "\n",
    "u = Function(V)\n",
    "bc = DirichletBC(V, 0, 'on_boundary')\n",
    "solve(lhs(a) == rhs(a), u, bcs=bc)\n",
    "\n",
    "c = trisurf(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14c95163-435d-4588-91a2-e2007ab4a31e",
   "metadata": {},
   "source": [
    "## Extrude for pml"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cec4de8-faa0-4b4e-93a6-357c3e43e5d0",
   "metadata": {},
   "source": [
    "### 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ad482da-c8e8-466c-95a5-a23eac3285e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import gmsh\n",
    "import numpy as np\n",
    "\n",
    "gmsh.initialize()\n",
    "sphere = gmsh.model.occ.add_sphere(0, 0, 0, 0.5, angle1=0)\n",
    "gmsh.model.occ.synchronize()\n",
    "bdy = gmsh.model.get_boundary([[3, sphere]])\n",
    "plane = []\n",
    "for dim, tag in bdy:\n",
    "    name = gmsh.model.get_type(dim, abs(tag))\n",
    "    if name == 'Sphere':\n",
    "        inner_bdy = tag\n",
    "    elif name == 'Plane':\n",
    "        plane.append(abs(tag))\n",
    "        \n",
    "# heights = [0.1 for _ in range(8)]\n",
    "# change recombine to True to get mix-cell mesh\n",
    "top = gmsh.model.geo.extrude_boundary_layer([[2, inner_bdy]], numElements=[2], heights=[0.1], recombine=False)\n",
    "\n",
    "gmsh.model.geo.synchronize()\n",
    "\n",
    "gmsh.model.add_physical_group(2, plane, tag=1)\n",
    "gmsh.model.set_physical_name(2, 1, \"plane\")\n",
    "gmsh.model.add_physical_group(2, [inner_bdy], tag=2)\n",
    "gmsh.model.set_physical_name(2, 2, \"inner bdy\")\n",
    "\n",
    "gmsh.model.add_physical_group(2, [top[3][1]], tag=3)\n",
    "gmsh.model.set_physical_name(2, 3, \"plane2\")\n",
    "gmsh.model.add_physical_group(2, [top[0][1]], tag=4)\n",
    "gmsh.model.set_physical_name(2, 4, \"outer bdy\")\n",
    "\n",
    "gmsh.model.add_physical_group(3, [sphere], tag=1)\n",
    "gmsh.model.set_physical_name(3, 1, \"Domain\")\n",
    "gmsh.model.add_physical_group(3, [top[1][1]], tag=2)\n",
    "gmsh.model.set_physical_name(3, 2, \"PML\")\n",
    "\n",
    "gmsh.model.geo.synchronize()\n",
    "# for dim, tag in top:\n",
    "#     name = gmsh.model.get_type(dim, abs(tag))\n",
    "#     print(name, dim, tag)\n",
    "    \n",
    "# gmsh.fltk.run()\n",
    "\n",
    "gmsh.model.mesh.generate()\n",
    "\n",
    "# gmsh.write('gmsh/sphere_pml.msh')\n",
    "\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a43218e-4183-430c-8267-787ff7ac9b6b",
   "metadata": {},
   "source": [
    "### 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b835811c-0f6a-43b4-9ce9-9ee25587b511",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import gmsh\n",
    "import numpy as np\n",
    "\n",
    "gmsh.initialize()\n",
    "disk = gmsh.model.occ.add_disk(0, 0, 0, 0.5, 0.5)\n",
    "\n",
    "gmsh.model.occ.synchronize()\n",
    "bdy = gmsh.model.get_boundary([[2, disk]])\n",
    "\n",
    "inner_bdy = bdy[0][1]\n",
    "\n",
    "       \n",
    "# heights = [0.1 for _ in range(8)]\n",
    "# change recombine to True to get mix-cell mesh\n",
    "top = gmsh.model.geo.extrude_boundary_layer([[1, inner_bdy]], numElements=[2], heights=[0.1], recombine=False)\n",
    "\n",
    "gmsh.model.geo.synchronize()\n",
    "\n",
    "gmsh.model.add_physical_group(1, [inner_bdy], tag=1)\n",
    "gmsh.model.set_physical_name(1, 1, \"inner bdy\")\n",
    "\n",
    "gmsh.model.add_physical_group(1, [top[0][1]], tag=2)\n",
    "gmsh.model.set_physical_name(1, 2, \"outer bdy\")\n",
    "\n",
    "gmsh.model.add_physical_group(2, [disk], tag=1)\n",
    "gmsh.model.set_physical_name(2, 1, \"Domain\")\n",
    "gmsh.model.add_physical_group(2, [top[1][1]], tag=2)\n",
    "gmsh.model.set_physical_name(2, 2, \"PML\")\n",
    "\n",
    "gmsh.model.geo.synchronize()\n",
    "# for dim, tag in top:\n",
    "#     name = gmsh.model.get_type(dim, abs(tag))\n",
    "#     print(name, dim, tag)\n",
    "    \n",
    "# gmsh.fltk.run()\n",
    "\n",
    "gmsh.model.mesh.generate()\n",
    "\n",
    "# gmsh.write('gmsh/disk_pml.msh')\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39a98150-c038-4cf6-8881-67e342d00757",
   "metadata": {},
   "source": [
    "#### Plot domain marker in space HDivT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349ab1e1-9a99-434a-8860-7b72bf9afdd7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from firedrake import *\n",
    "import numpy as np\n",
    "\n",
    "mesh = Mesh('gmsh/disk_pml.msh')\n",
    "V = FunctionSpace(mesh, 'HDivT', 1)\n",
    "\n",
    "\n",
    "marker = Function(V, name='f')\n",
    "par_loop(('{[i] : 0 <= i < f.dofs}', 'f[i, 0] = 1/2'), dx(1), {'f': (marker, INC)})\n",
    "\n",
    "index = marker.dat.data_with_halos < 3/4\n",
    "marker.dat.data_with_halos[index] = 0\n",
    "\n",
    "plex = mesh.topology_dm\n",
    "\n",
    "s, e = plex.getHeightStratum(1)\n",
    "\n",
    "coords = mesh.coordinates\n",
    "csec = coords.function_space().dm.getSection()\n",
    "sec = V.dm.getSection()\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=[8, 8])\n",
    "for i in range(s, e):\n",
    "    a, b = plex.getCone(i)\n",
    "\n",
    "    off_a = csec.getOffset(a)\n",
    "    off_b = csec.getOffset(b)\n",
    "    x1, y1 = coords.dat.data_ro_with_halos[off_a].real\n",
    "    x2, y2 = coords.dat.data_ro_with_halos[off_b].real\n",
    "    plt.plot([x1, x2], [y1, y2])\n",
    "    off_i = sec.getOffset(i)\n",
    "    v = marker.dat.data_ro_with_halos[off_i].real\n",
    "    plt.text((x1 + x2)/2, (y1 + y2)/2, round(v), ha='center', va='center')\n",
    "\n",
    "plt.axis('equal')\n",
    "\n",
    "rank, size = mesh.comm.rank, mesh.comm.size\n",
    "plt.savefig(f'figures/hdivt-marker-{size}-{rank}.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aab0bcf0-0d5a-4694-818c-459957176621",
   "metadata": {},
   "source": [
    "# Gmsh Simplex Ordering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "371d9828-75c9-4ab1-a814-db2add21d500",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def SN1(p):\n",
    "    return p + 1\n",
    "\n",
    "def SN2(p):\n",
    "    return SN1(p) * SN1(p + 1) // 2\n",
    "    \n",
    "def SN3(p): \n",
    "    return SN2(p) * SN1(p + 2) // 3\n",
    "    \n",
    "def SI1(p, i):\n",
    "    return i\n",
    "\n",
    "def SI2(p, i, j): \n",
    "    return i + (SN2(p) - SN2(p - j))\n",
    "\n",
    "def SI3(p, i, j, k):\n",
    "    return SI2(p - k, i, j) + SN3(p) - SN3(p - k)\n",
    "\n",
    "def SL1(p):\n",
    "    for i in range(1, p):\n",
    "        yield i\n",
    "\n",
    "def SL2(p):\n",
    "    for i in range(1, p - 1):\n",
    "        for j in range(1, p - i):\n",
    "            yield i, j\n",
    "\n",
    "def SL3(p):\n",
    "    for i in range(1, p - 2):\n",
    "        for j in range(1, p - i):\n",
    "            for k in range(1, p - i - j):\n",
    "                yield i, j, k\n",
    "\n",
    "def GmshLexOrder_SEG(p, node=0):\n",
    "    index = lambda i: SI1(p, i)\n",
    "\n",
    "    lex = - np.ones(SN1(p))\n",
    "    \n",
    "    if p == 0:\n",
    "        lex[0] = node; node += 1\n",
    "        return lex, node\n",
    "    \n",
    "    lex[index(0)] = node; node += 1\n",
    "    lex[index(p)] = node; node += 1\n",
    "    if p == 1:\n",
    "        return lex, node\n",
    "    \n",
    "    for i in SL1(p):\n",
    "        lex[index(i)] = node; node += 1\n",
    "        \n",
    "    return lex, node\n",
    "\n",
    "def GmshLexOrder_TRI(p, node=0):\n",
    "    index = lambda i, j: SI2(p, i, j)\n",
    "    \n",
    "    lex = - np.ones(SN2(p))\n",
    "    \n",
    "    if p == 0:\n",
    "        lex[0] = node; node += 1\n",
    "        return lex, node\n",
    "    \n",
    "    lex[index(0, 0)] = node; node += 1\n",
    "    lex[index(p, 0)] = node; node += 1\n",
    "    lex[index(0, p)] = node; node += 1\n",
    "    \n",
    "    if p == 1:\n",
    "        return lex, node\n",
    "\n",
    "    for i in SL1(p):\n",
    "        lex[index(i, 0)]     = node; node += 1\n",
    "    for j in SL1(p):\n",
    "        lex[index(p - j, j)] = node; node += 1\n",
    "    for j in SL1(p):\n",
    "        lex[index(0, p - j)] = node; node += 1\n",
    "\n",
    "    if p == 2:\n",
    "        return lex, node\n",
    "    \n",
    "    sub, node = GmshLexOrder_TRI(p - 3, node);\n",
    "    for _, (j, i) in enumerate(SL2(p)):\n",
    "        lex[index(i, j)] = sub[_];\n",
    "\n",
    "    return lex, node\n",
    "\n",
    "def GmshLexOrder_TET(p, node=0):\n",
    "    index = lambda i, j, k: SI3(p, i, j, k)\n",
    "    lex = - np.ones(SN3(p))\n",
    "    \n",
    "    if p == 0:\n",
    "        lex[0] = node; node += 1\n",
    "        return lex, node\n",
    "    lex[index(0, 0, 0)] = node; node += 1\n",
    "    lex[index(p, 0, 0)] = node; node += 1\n",
    "    lex[index(0, p, 0)] = node; node += 1\n",
    "    lex[index(0, 0, p)] = node; node += 1\n",
    "    \n",
    "    if p == 1:\n",
    "        return lex, node\n",
    "    \n",
    "    # internal edge nodes \n",
    "    for i in SL1(p): lex[index(    i,     0,     0)] = node; node += 1\n",
    "    for j in SL1(p): lex[index(p - j,     j,     0)] = node; node += 1\n",
    "    for j in SL1(p): lex[index(    0, p - j,     0)] = node; node += 1\n",
    "    for k in SL1(p): lex[index(    0,     0, p - k)] = node; node += 1\n",
    "    for j in SL1(p): lex[index(    0,     j, p - j)] = node; node += 1\n",
    "    for i in SL1(p): lex[index(    i,     0, p - i)] = node; node += 1\n",
    "    \n",
    "    if p == 2:\n",
    "        return lex, node\n",
    "    \n",
    "    # /* internal face nodes */\n",
    "    sub, node = GmshLexOrder_TRI(p - 3, node)\n",
    "    for _, (i, j) in enumerate(SL2(p)): \n",
    "        lex[index(i, j, 0)] = sub[_]\n",
    "        \n",
    "    sub, node = GmshLexOrder_TRI(p - 3, node);\n",
    "    for _, (k, i) in enumerate(SL2(p)): \n",
    "        lex[index(i, 0, k)] = sub[_]\n",
    "        \n",
    "    sub, node = GmshLexOrder_TRI(p - 3, node);\n",
    "    for _, (j, k) in enumerate(SL2(p)): \n",
    "        lex[index(0, j, k)] = sub[_]\n",
    "        \n",
    "    sub, node = GmshLexOrder_TRI(p - 3, node);\n",
    "    for _, (j, i) in enumerate(SL2(p)): \n",
    "        lex[index(i, j, p - i - j)] = sub[_]\n",
    "        \n",
    "    if p == 3: \n",
    "        return lex, node\n",
    "\n",
    "    # internal cell nodes */\n",
    "    sub, node = GmshLexOrder_TET(p - 4, node);\n",
    "    for _, (k, j, i) in enumerate(SL3(p)):\n",
    "        lex[index(i, j, k)] = sub[_];\n",
    "\n",
    "    return lex, node"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd37b6b2-4031-4371-8e31-d6ac8c3cdc20",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Test high order\n",
    "\n",
    "TODO: the code do not work well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae140879-12f5-41a6-b544-3c561e830331",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from firedrake import *\n",
    "import firedrake as fd\n",
    "import ufl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "try:\n",
    "    from petsctools.options import OptionsManager\n",
    "except ImportError:\n",
    "    from firedrake.petsc import OptionsManager"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8c1b978-51b0-4c51-828f-f64e2ed58d78",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def getCoordinateFESpaceOrder(dm):\n",
    "    cdm = dm.getCoordinateDM()\n",
    "    kls, _ = cdm.getField(0)\n",
    "    if kls.getClassName() == 'PetscFE':\n",
    "        p = int(kls.getName()[1:])\n",
    "    else:\n",
    "        p = 1\n",
    "    return p\n",
    "\n",
    "def callback(mesh):\n",
    "\n",
    "    \"\"\"Finish initialisation.\"\"\"\n",
    "    del mesh._callback\n",
    "\n",
    "    mesh.topology.init()\n",
    "\n",
    "    coordinates_fs = fd.functionspace.FunctionSpace(mesh.topology, mesh.ufl_coordinate_element())\n",
    "    sec = coordinates_fs.dm.getDefaultSection()\n",
    "\n",
    "    dm = mesh.topology.topology_dm\n",
    "    dim = dm.getCoordinateDim()\n",
    "    dm_sec = dm.getCoordinateSection()\n",
    "    dm_coords = dm.getCoordinatesLocal().array_r\n",
    "    \n",
    "    coordinates_data = np.empty_like(dm_coords).reshape([-1, dim])\n",
    "    s, e = sec.getChart()\n",
    "    for i in range(s, e):\n",
    "        ndof = sec.getDof(i)\n",
    "        if ndof > 0:\n",
    "            offset = sec.getOffset(i)\n",
    "            dm_offset = dm_sec.getOffset(i)\n",
    "            coordinates_data[offset, :] = dm_coords[dm_offset:(dm_offset+dim)]\n",
    "\n",
    "    # Finish the initialisation of mesh topology\n",
    "    coordinates = fd.function.CoordinatelessFunction(coordinates_fs, val=coordinates_data, name=mesh.name + \"_coordinates\")\n",
    "\n",
    "    mesh.__init__(coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2554b347-1302-4019-80ad-fa7c6752c3b7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# https://github.com/firedrakeproject/firedrake/blob/ec0329f092b431e8e4c8bd7e41f6667234c9caa3/firedrake/mesh.py#L2827\n",
    "def make_mesh_from_mesh_topology(topology, name):\n",
    "    import finat.ufl\n",
    "    # Construct coordinate element\n",
    "    # TODO: meshfile might indicates higher-order coordinate element\n",
    "    cell = topology.ufl_cell()\n",
    "    geometric_dim = topology.topology_dm.getCoordinateDim()\n",
    "    cell = cell.reconstruct(geometric_dimension=geometric_dim)\n",
    "    topology_dim = topology.topology_dm.getDimension()\n",
    "    \n",
    "    p = getCoordinateFESpaceOrder(topology.topology_dm)\n",
    "    # TODO: we only process Lagrange Now\n",
    "    element = finat.ufl.VectorElement(\"Lagrange\", cell, p)\n",
    "    mesh = fd.mesh.MeshGeometry.__new__(fd.mesh.MeshGeometry, element, topology.comm)\n",
    "    mesh._init_topology(topology)\n",
    "    mesh.name = name\n",
    "    if p > 1:\n",
    "        mesh._callback = callback\n",
    "    return mesh\n",
    "\n",
    "fd.mesh.make_mesh_from_mesh_topology = make_mesh_from_mesh_topology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbe2c2c8-d56b-4769-ab92-3881576ff83c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "project_parameter = {\n",
    "    \"dm_plex_gmsh_project\": True,\n",
    "    \"dm_plex_gmsh_project_\": {\n",
    "        # \"fe_view\": \"ascii\",\n",
    "        \"petscdualspace_lagrange_continuity\": True,\n",
    "    }\n",
    "}\n",
    "\n",
    "om = OptionsManager(project_parameter, options_prefix=\"\")\n",
    "with om.inserted_options():\n",
    "    mesh = Mesh('gmsh/cube_2rd.msh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa39ce00-da0b-4f92-bc89-df34a7b8e0f8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "c = triplot(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16510ec2-7af3-476b-8f9e-224efe3d77fd",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Show the issue of gmsh read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a5c0943-00c1-4e0d-b5a3-a76c7386fabe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from firedrake.petsc import PETSc\n",
    "import numpy as np\n",
    "\n",
    "try:\n",
    "    from petsctools.options import OptionsManager\n",
    "except ImportError:\n",
    "    from firedrake.petsc import OptionsManager\n",
    "\n",
    "def show_plex_coordinates_continue(plex):\n",
    "    dim = plex.getCoordinateDim()\n",
    "    dm_sec = plex.getCoordinateSection()\n",
    "    dm_coords = plex.getCoordinatesLocal().array_r.real.reshape([-1, dim])\n",
    "    cs, ce = plex.getHeightStratum(0)\n",
    "\n",
    "    print(f\"plex {plex.getName()}\")\n",
    "    for i in range(cs, ce):\n",
    "        fs = plex.getCone(i)\n",
    "        print(\"  cell:\", i, tuple(fs))\n",
    "        cl, o = plex.getTransitiveClosure(i)\n",
    "        print(\"    coordinates:\")\n",
    "        for k, p in enumerate(cl[6:]):\n",
    "            offset = dm_sec.getOffset(p)\n",
    "            print(\"      \", p, \"[%.2f, %.2f, %.2f]\"%tuple(dm_coords[offset//dim, :]))\n",
    "            \n",
    "def show_plex_coordinates(plex):\n",
    "    dim = plex.getCoordinateDim()\n",
    "    dm_sec = plex.getCoordinateSection()\n",
    "    dm_coords = plex.getCoordinatesLocal().array_r.real\n",
    "    cs, ce = plex.getHeightStratum(0)\n",
    "\n",
    "    print(f\"plex {plex.getName()}\")\n",
    "    for i in range(cs, ce):\n",
    "        offset = dm_sec.getOffset(i)\n",
    "        dof = dm_sec.getDof(i)\n",
    "        fs = plex.getCone(i)\n",
    "\n",
    "        print(\"  cell:\", i, tuple(fs))\n",
    "        print(\"    coordinates:\")\n",
    "        coords = dm_coords[offset:offset+dof].reshape([-1, dim])\n",
    "        for k in coords:\n",
    "            print(\"      \", \"[%.2f, %.2f, %.2f]\"%tuple(k))\n",
    "\n",
    "plex = PETSc.DMPlex().createFromFile('gmsh/cube_2rd.msh', 'gmsh_2rd')\n",
    "show_plex_coordinates(plex)\n",
    "\n",
    "project_parameter = {\n",
    "    \"dm_plex_gmsh_project\": True,\n",
    "    \"dm_plex_gmsh_project_\": {\n",
    "        # \"fe_view\": \"ascii\",\n",
    "        \"petscdualspace_lagrange_continuity\": True,\n",
    "    }\n",
    "}\n",
    "\n",
    "om = OptionsManager(project_parameter, options_prefix=\"\")\n",
    "with om.inserted_options():\n",
    "    plex_proj = PETSc.DMPlex().createFromFile('gmsh/cube_2rd.msh', 'gmsh_2rd_with_projection')\n",
    "show_plex_coordinates_continue(plex_proj)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}