Improvements to skfem.io.meshio

skfem/io/meshio.py

@@ -19,45 +19,72 @@
     'quad9': skfem.MeshQuad2,
 }
 
+BOUNDARY_TYPE_MAPPING = {
+    'line': 'vertex',
+    'triangle': 'line',
+    'quad': 'line',
+    'tetra': 'triangle',
+    'hexahedron': 'quad',
+    'tetra10': 'triangle',  # TODO support quadratic facets
+    'triangle6': 'line',  # TODO
+    'quad9': 'line',  # TODO
+}
+
 TYPE_MESH_MAPPING = {MESH_TYPE_MAPPING[k]: k
                      for k in dict(reversed(list(MESH_TYPE_MAPPING.items())))}
 
 
-def from_meshio(m, force_mesh_type=None):
+def from_meshio(m,
+                int_data_to_sets=False):
     """Convert meshio mesh into :class:`skfem.mesh.Mesh`.
 
     Parameters
     ----------
     m
         The mesh from meshio.
-    force_mesh_type
-        An optional string forcing the mesh type if automatic detection
-        fails. See :data:`skfem.io.meshio.MESH_TYPE_MAPPING` for possible
-        values.
+    int_data_to_sets
+        We correctly read the so-called "cell sets" from ``meshio``.  However,
+        many mesh formats do not support "cell sets" natively and, instead, use
+        cellwise integer data to distinguish between different subdomains and
+        boundaries.  If ``True``, call ``meshio.Mesh.int_data_to_sets`` to
+        convert between the representations before attempting to read tags from
+        ``meshio``.
 
     Returns
     -------
     A :class:`~skfem.mesh.Mesh` object.
 
     """
-
     cells = m.cells_dict
 
-    if force_mesh_type is None:
-        meshio_type = None
 
-        for k, v in MESH_TYPE_MAPPING.items():
-            # find first if match
-            if k in cells:
-                meshio_type, mesh_type = k, v
+    meshio_type = None
+
+    # detect 3D
+    for k in cells:
+        if k in {'tetra', 'hexahedron', 'tetra10'}:
+            meshio_type = k
+            break
+
+    if meshio_type is None:
+        # detect 2D
+        for k in cells:
+            if k in {'triangle', 'quad', 'triangle6', 'quad9'}:
+                meshio_type = k
                 break
 
-        if meshio_type is None:
-            raise NotImplementedError("Mesh type(s) not supported "
-                                      "in import: {}.".format(cells.keys()))
-    else:
-        meshio_type, mesh_type = (force_mesh_type,
-                                  MESH_TYPE_MAPPING[force_mesh_type])
+    if meshio_type is None:
+        # detect 1D
+        for k in cells:
+            if k == 'line':
+                meshio_type = k
+                break
+
+    if meshio_type is None:
+        raise NotImplementedError("Mesh type(s) not supported "
+                                  "in import: {}.".format(cells.keys()))
+
+    mesh_type = MESH_TYPE_MAPPING[meshio_type]
 
     # create p and t
     p = np.ascontiguousarray(mesh_type.strip_extra_coordinates(m.points).T)
@@ -70,21 +97,17 @@ def from_meshio(m, force_mesh_type=None):
     mtmp = mesh_type(p, t)
 
     try:
-        # element to boundary element type mapping
-        bnd_type = {
-            'line': 'vertex',
-            'triangle': 'line',
-            'quad': 'line',
-            'tetra': 'triangle',
-            'hexahedron': 'quad',
-        }[meshio_type]
+        bnd_type = BOUNDARY_TYPE_MAPPING[meshio_type]
 
         def find_tagname(tag):
             for key in m.field_data:
                 if m.field_data[key][0] == tag:
                     return key
             return None
 
+        if int_data_to_sets:
+            m.int_data_to_sets()
+
         if m.cell_sets:  # MSH 4.1
             subdomains = {k: v[meshio_type]
                           for k, v in m.cell_sets_dict.items()
@@ -137,19 +160,71 @@ def find_tagname(tag):
     return mtmp
 
 
-def from_file(filename):
-    return from_meshio(meshio.read(filename))
+def from_file(filename, **kwargs):
+    return from_meshio(meshio.read(filename), **kwargs)
 
 
-def to_meshio(mesh, point_data=None):
+def to_meshio(mesh,
+              point_data=None,
+              cell_data=None,
+              write_boundaries=False,
+              sets_to_int_data=False):
 
     t = mesh.dofs.element_dofs.copy()
     if isinstance(mesh, skfem.MeshHex):
         t = t[[0, 3, 6, 2, 1, 5, 7, 4]]
 
-    cells = {TYPE_MESH_MAPPING[type(mesh)]: t.T}
-    return meshio.Mesh(mesh.p.T, cells, point_data)
-
-
-def to_file(mesh, filename, point_data=None, **kwargs):
-    meshio.write(filename, to_meshio(mesh, point_data), **kwargs)
+    mtype = TYPE_MESH_MAPPING[type(mesh)]
+    bmtype = BOUNDARY_TYPE_MAPPING[mtype]
+    cells = {mtype: t.T}
+
+    if write_boundaries:
+        bfacets = mesh.boundary_facets()
+        cells[bmtype] = mesh.facets[:, bfacets].T
+
+    # write named subsets
+    cell_sets = {}
+    if mesh.subdomains is not None:
+        for key in mesh.subdomains:
+            cell_sets[key] = [
+                mesh.subdomains[key],
+                *((len(cells) - 1) * (np.array([], dtype=np.int64),)),
+            ]
+
+    if write_boundaries and mesh.boundaries is not None:
+        ix_map = np.zeros(mesh.facets.shape[1], dtype=np.int64) - 1
+        ix_map[bfacets] = np.arange(len(bfacets), dtype=np.int64)
+        for key in mesh.boundaries:
+            cell_sets[key] = [
+                np.array([], dtype=np.int64),
+                ix_map[mesh.boundaries[key]],
+            ]
+
+    mio = meshio.Mesh(
+        mesh.p.T,
+        cells,
+        point_data=point_data,
+        cell_data=cell_data,
+        cell_sets=cell_sets,
+    )
+
+    if sets_to_int_data:
+        mio.sets_to_int_data()
+
+    return mio
+
+
+def to_file(mesh,
+            filename,
+            point_data=None,
+            cell_data=None,
+            write_boundaries=False,
+            sets_to_int_data=False,
+            **kwargs):
+    meshio.write(filename,
+                 to_meshio(mesh,
+                           point_data,
+                           cell_data,
+                           write_boundaries,
+                           sets_to_int_data),
+                 **kwargs)

skfem/mesh/mesh.py

@@ -407,6 +407,9 @@ def __str__(self):
     def save(self,
              filename: str,
              point_data: Optional[Dict[str, ndarray]] = None,
+             cell_data: Optional[Dict[str, ndarray]] = None,
+             write_boundaries: bool = False,
+             sets_to_int_data: bool = False,
              **kwargs) -> None:
         """Export the mesh and fields using meshio.
 
@@ -420,10 +423,16 @@ def save(self,
 
         """
         from skfem.io.meshio import to_file
-        return to_file(self, filename, point_data, **kwargs)
+        return to_file(self,
+                       filename,
+                       point_data,
+                       cell_data,
+                       write_boundaries,
+                       sets_to_int_data,
+                       **kwargs)
 
     @classmethod
-    def load(cls, filename: str):
+    def load(cls, filename: str, **kwargs):
         """Load a mesh using meshio.
 
         Parameters
@@ -433,7 +442,7 @@ def load(cls, filename: str):
 
         """
         from skfem.io.meshio import from_file
-        return from_file(filename)
+        return from_file(filename, **kwargs)
 
     @classmethod
     def from_dict(cls, data):

tests/test_mesh.py

@@ -281,6 +281,7 @@ def test_finder_simplex(m, seed):
         MeshTri2(),
         MeshQuad2(),
         MeshTet2(),
+        MeshLine(),
     ]
 )
 def test_meshio_cycle(m):
@@ -290,26 +291,98 @@ def test_meshio_cycle(m):
     assert_array_equal(M.t, m.t)
 
 
+_test_lambda = {
+    'left': lambda x: x[0] < 0.5,
+    'right': lambda x: x[0] > 0.5,
+}
+
+
 @pytest.mark.parametrize(
     "m",
     [
         MeshTri(),
         MeshQuad(),
+        # MeshHex(),  # TODO facet order changes?
         MeshTet(),
+    ]
+)
+def test_meshio_cycle_boundaries(m):
+
+    m = m.with_boundaries(_test_lambda)
+    M = from_meshio(to_meshio(m, write_boundaries=True))
+    assert_array_equal(M.p, m.p)
+    assert_array_equal(M.t, m.t)
+    for key in m.boundaries:
+        assert_array_equal(M.boundaries[key], m.boundaries[key])
+
+
+@pytest.mark.parametrize(
+    "m",
+    [
+        MeshTri(),
+        MeshQuad(),
         MeshHex(),
+        MeshTet(),
     ]
 )
-def test_saveload_cycle(m):
+def test_meshio_cycle_subdomains(m):
+
+    m = m.refined().with_subdomains(_test_lambda)
+    M = from_meshio(to_meshio(m))
+    assert_array_equal(M.p, m.p)
+    assert_array_equal(M.t, m.t)
+    for key in m.subdomains:
+        assert_array_equal(M.subdomains[key], m.subdomains[key])
+
+
+@pytest.mark.parametrize(
+    "m",
+    [
+        MeshTri(),
+        MeshQuad(),
+        MeshTet(),
+        MeshHex(),
+    ]
+)
+def test_saveload_cycle_vtk(m):
 
     from tempfile import NamedTemporaryFile
     m = m.refined(2)
-    f = NamedTemporaryFile(delete=False)
-    m.save(f.name + ".vtk")
-    with pytest.warns(UserWarning):
-       m2 = Mesh.load(f.name + ".vtk")
+    with NamedTemporaryFile() as f:
+        m.save(f.name + ".vtk")
+        with pytest.warns(UserWarning):
+            m2 = Mesh.load(f.name + ".vtk")
+
+    assert_array_equal(m.p, m2.p)
+    assert_array_equal(m.t, m2.t)
+
+
+@pytest.mark.parametrize(
+    "m",
+    [
+        MeshTri(),
+        MeshQuad(),
+        # MeshHex(),  # TODO facet order changes?
+        MeshTet(),
+    ]
+)
+def test_saveload_cycle_tags(m):
+
+    m = (m
+         .refined()
+         .with_subdomains(_test_lambda)
+         .with_boundaries({'test': lambda x: x[0] == 0}))
+    from tempfile import NamedTemporaryFile
+    with NamedTemporaryFile() as f:
+        m.save(f.name + ".inp", write_boundaries=True)
+        m2 = Mesh.load(f.name + ".inp")
 
     assert_array_equal(m.p, m2.p)
     assert_array_equal(m.t, m2.t)
+    for key in m.subdomains:
+        assert_array_equal(m2.subdomains[key], m.subdomains[key])
+    for key in m.boundaries:
+        assert_array_equal(m2.boundaries[key], m.boundaries[key])
 
 
 if __name__ == '__main__':