Merge pull request #365 from nschloe/generate-to-mesh

Generate to mesh

README.md

@@ -37,7 +37,7 @@ Codes:
 ```python
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     geom.add_polygon(
         [
             [0.0, 0.0],
@@ -47,22 +47,22 @@ with pygmsh.built_in.Geometry() as geom:
         ],
         mesh_size=0.1,
     )
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 
 # mesh.points, mesh.cells, ...
 # mesh.write("out.vtk")
 ```
 ```python
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 ```python
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     lcar = 0.1
     p1 = geom.add_point([0.0, 0.0, 0.0], lcar)
     p2 = geom.add_point([1.0, 0.0, 0.0], lcar)
@@ -77,7 +77,7 @@ with pygmsh.built_in.Geometry() as geom:
     ll = geom.add_curve_loop([s1, s2])
     pl = geom.add_plane_surface(ll)
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 
 The return value is always a [meshio](https://pypi.org/project/meshio) mesh, so to store
@@ -98,7 +98,7 @@ The output file can be visualized with various tools, e.g.,
 ```python
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     poly = geom.add_polygon(
         [
             [0.0, 0.0],
@@ -109,13 +109,13 @@ with pygmsh.built_in.Geometry() as geom:
         mesh_size=0.1,
     )
     geom.extrude(poly, [0.0, 0.3, 1.0], num_layers=5)
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 ```python
 from math import pi
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     poly = geom.add_polygon(
         [
             [0.0, 0.2, 0.0],
@@ -125,13 +125,13 @@ with pygmsh.built_in.Geometry() as geom:
         mesh_size=0.1,
     )
     geom.revolve(poly, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0], 0.8 * pi)
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 ```python
 from math import pi
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     poly = geom.add_polygon(
         [
             [+0.0, +0.5],
@@ -154,21 +154,21 @@ with pygmsh.built_in.Geometry() as geom:
         angle=pi / 3,
     )
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 
 #### OpenCASCADE
 <img src="https://nschloe.github.io/pygmsh/intersection.png" width="100%"> | <img src="https://nschloe.github.io/pygmsh/ellipsoid-holes.png" width="100%"> | <img src="https://nschloe.github.io/pygmsh/puzzle.png" width="100%">
 :------------------:|:-------------:|:--------:|
  |    |   |
 
-As of version 3.0, Gmsh supports OpenCASCADE, allowing for a CAD-style geometry
+As of version 3.0, Gmsh supports OpenCASCADE (`occ`), allowing for a CAD-style geometry
 specification.
 ```python
 from math import pi, cos
 import pygmsh
 
-with pygmsh.opencascade.Geometry() as geom:
+with pygmsh.occ.Geometry() as geom:
     geom.characteristic_length_max = 0.1
     r = 0.5
     disks = [
@@ -178,13 +178,13 @@ with pygmsh.opencascade.Geometry() as geom:
     ]
     geom.boolean_intersection(disks)
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 ```python
 # ellpsoid with holes
 import pygmsh
 
-with pygmsh.opencascade.Geometry() as geom:
+with pygmsh.occ.Geometry() as geom:
     geom.characteristic_length_max = 0.1
     ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [1.0, 0.7, 0.5])
 
@@ -195,13 +195,13 @@ with pygmsh.opencascade.Geometry() as geom:
     ]
     geom.boolean_difference(ellipsoid, geom.boolean_union(cylinders))
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 ```python
 # puzzle piece
 import pygmsh
 
-with pygmsh.opencascade.Geometry() as geom:
+with pygmsh.occ.Geometry() as geom:
     geom.characteristic_length_min = 0.1
     geom.characteristic_length_max = 0.1
 
@@ -218,7 +218,7 @@ with pygmsh.opencascade.Geometry() as geom:
 
     geom.extrude(flat, [0, 0, 0.3])
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 
 
@@ -231,7 +231,7 @@ with pygmsh.opencascade.Geometry() as geom:
 # boundary refinement
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     poly = geom.add_polygon(
         [
             [0.0, 0.0, 0.0],
@@ -259,14 +259,14 @@ with pygmsh.built_in.Geometry() as geom:
     )
     geom.set_background_mesh([field0, field1], operator="Min")
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 <!--exdown-skip-->
 ```python
 # mesh refinement with callback
 import pygmsh
 
-with pygmsh.built_in.Geometry() as geom:
+with pygmsh.geo.Geometry() as geom:
     geom.add_polygon(
         [
             [-1.0, -1.0],
@@ -279,7 +279,7 @@ with pygmsh.built_in.Geometry() as geom:
         lambda dim, tag, x, y, z: 6.0e-2 + 2.0e-1 * (x ** 2 + y ** 2)
     )
 
-    mesh = pygmsh.generate_mesh(geom)
+    mesh = geom.generate_mesh()
 ```
 <!--exdown-skip-->
 ```python
@@ -288,7 +288,7 @@ from math import sqrt
 import pygmsh
 
 
-with pygmsh.opencascade.Geometry() as geom:
+with pygmsh.occ.Geometry() as geom:
     geom.add_ball([0.0, 0.0, 0.0], 1.0)
 
     geom.set_mesh_size_callback(

pygmsh/__init__.py

@@ -1,11 +1,10 @@
-from . import built_in, opencascade
+from . import geo, occ
 from .__about__ import __gmsh_version__, __version__
-from .helpers import generate_mesh, orient_lines, rotation_matrix
+from .helpers import orient_lines, rotation_matrix
 
 __all__ = [
-    "built_in",
-    "opencascade",
-    "generate_mesh",
+    "geo",
+    "occ",
     "rotation_matrix",
     "orient_lines",
     "__version__",

pygmsh/common/geometry.py

@@ -1,6 +1,8 @@
 import math
 
 import gmsh
+import meshio
+import numpy
 
 from .bspline import BSpline
 from .circle_arc import CircleArc
@@ -19,7 +21,7 @@
 
 class CommonGeometry:
     """Geometry base class containing all methods that can be shared between built-in
-    and opencascade.
+    and occ.
     """
 
     def __init__(self, env):
@@ -301,3 +303,137 @@ def in_volume(self, input_entity, volume):
 
     def set_mesh_size_callback(self, fun):
         gmsh.model.mesh.setSizeCallback(fun)
+
+    def generate_mesh(
+        self,
+        dim=3,
+        order=None,
+        prune_vertices=True,
+        prune_z_0=False,
+        remove_lower_dim_cells=False,
+    ):
+        """Return a meshio.Mesh, storing the mesh points, cells, and data, generated by
+        Gmsh from the `self`, written to a temporary file, and reread by `meshio`.
+
+        Gmsh's native "msh" format is ill-suited to fast I/O.  This can greatly reduce
+        the performance of pygmsh.  As alternatives, try `mesh_file_type=`:
+
+        - "vtk"`, though Gmsh doesn't write the physical tags to VTK
+          <https://gitlab.onelab.info/gmsh/gmsh/issues/389> or
+
+        - `"mesh"`, though this only supports a few basic elements - "line", "triangle",
+          "quad", "tetra", "hexahedron" - and doesn't preserve the `$PhysicalNames`,
+          just the `int` tags.
+        """
+        self.synchronize()
+
+        for item in self._AFTER_SYNC_QUEUE:
+            item.exec()
+
+        for item, host in self._EMBED_QUEUE:
+            gmsh.model.mesh.embed(item.dimension, [item._ID], host.dimension, host._ID)
+
+        # set compound entities after sync
+        for c in self._COMPOUND_ENTITIES:
+            gmsh.model.mesh.setCompound(*c)
+
+        for s in self._RECOMBINE_ENTITIES:
+            gmsh.model.mesh.setRecombine(*s)
+
+        for t in self._TRANSFINITE_CURVE_QUEUE:
+            gmsh.model.geo.mesh.setTransfiniteCurve(*t)
+
+        for t in self._TRANSFINITE_SURFACE_QUEUE:
+            gmsh.model.geo.mesh.setTransfiniteSurface(*t)
+
+        for item, size in self._SIZE_QUEUE:
+            gmsh.model.mesh.setSize(
+                gmsh.model.getBoundary(item.dim_tags, False, False, True), size
+            )
+
+        if order is not None:
+            gmsh.model.mesh.setOrder(order)
+
+        gmsh.model.mesh.generate(dim)
+
+        # extract point coords
+        idx, points, _ = gmsh.model.mesh.getNodes()
+        points = points.reshape(-1, 3)
+        idx -= 1
+        srt = numpy.argsort(idx)
+        assert numpy.all(idx[srt] == numpy.arange(len(idx)))
+        points = points[srt]
+        if prune_z_0 and numpy.all(numpy.abs(points[:, 2]) < 1.0e-13):
+            points = points[:, :2]
+
+        # extract cells
+        elem_types, elem_tags, node_tags = gmsh.model.mesh.getElements()
+        cells = []
+        for elem_type, node_tags in zip(elem_types, node_tags):
+            # `elementName', `dim', `order', `numNodes', `localNodeCoord', `numPrimaryNodes'
+            num_nodes_per_cell = gmsh.model.mesh.getElementProperties(elem_type)[3]
+            meshio.gmsh.gmsh_to_meshio_type
+            cells.append(
+                meshio.CellBlock(
+                    meshio.gmsh.gmsh_to_meshio_type[elem_type],
+                    node_tags.reshape(-1, num_nodes_per_cell) - 1,
+                )
+            )
+
+        # print("a", gmsh.model.getEntities())
+        # grps = gmsh.model.getPhysicalGroups()
+        # print("a", grps)
+        # for dim, tag in grps:
+        #     print("a", gmsh.model.getPhysicalName(dim, tag))
+        #     ent = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
+        #     print("a", ent)
+        #     assert len(ent) == 1
+        #     print("a", gmsh.model.mesh.getElements(dim, ent[0]))
+
+        # make meshio mesh
+        mesh = meshio.Mesh(points, cells)
+
+        if remove_lower_dim_cells:
+            # Only keep the cells of highest topological dimension; discard faces and such.
+            cells_2d = {"triangle", "quad"}
+            cells_3d = {
+                "tetra",
+                "hexahedron",
+                "wedge",
+                "pyramid",
+                "penta_prism",
+                "hexa_prism",
+            }
+            if any(c.type in cells_3d for c in mesh.cells):
+                keep_types = cells_3d
+            elif any(c.type in cells_2d for c in mesh.cells):
+                keep_types = cells_2d
+            else:
+                keep_types = set(cell_type for cell_type, _ in mesh.cells)
+
+            for name, val in mesh.cell_data.items():
+                mesh.cell_data[name] = [
+                    d for d, c in zip(val, mesh.cells) if c[0] in keep_types
+                ]
+            mesh.cells = [c for c in mesh.cells if c[0] in keep_types]
+
+        if prune_vertices:
+            # Make sure to include only those vertices which belong to a cell.
+            ncells = numpy.concatenate([numpy.concatenate(c) for _, c in mesh.cells])
+            uvertices, uidx = numpy.unique(ncells, return_inverse=True)
+
+            k = 0
+            cells = []
+            for key, cellblock in mesh.cells:
+                n = numpy.prod(cellblock.shape)
+                cells.append(
+                    meshio.CellBlock(key, uidx[k : k + n].reshape(cellblock.shape))
+                )
+                k += n
+            mesh.cells = cells
+
+            mesh.points = mesh.points[uvertices]
+            for key in mesh.point_data:
+                mesh.point_data[key] = mesh.point_data[key][uvertices]
+
+        return mesh