Test Meshes with Multiple Groups

At the moment this just adds a test for a mesh with 2 groups. Basic operations work just fine, but it currently fails in make_opposite_face_connection, so some fixes are expected.

meshmode/mesh/processing.py

@@ -37,6 +37,7 @@ __doc__ = """
 .. autofunction:: perform_flips
 .. autofunction:: find_bounding_box
 .. autofunction:: merge_disjoint_meshes
+.. autofunction:: split_mesh_groups
 .. autofunction:: map_mesh
 .. autofunction:: affine_map
 """
@@ -572,9 +573,19 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
 # {{{ split meshes
 
 
-def split_mesh_groups(mesh, element_flags):
-    """Split all the groups in *mesh* in two according to the values of
-    *element_flags*.
+def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
+    """Split all the groups in *mesh* in according to the values of
+    *element_flags*. The element flags are expected to be integers
+    defining, for each group, how the elements are to be split into
+    subgroups. For example, a single-group mesh with flags
+
+        .. code::
+
+            element_flags = [0, 0, 0, 42, 42, 42, 0, 0, 0, 41, 41, 41]
+
+    will create three subgroups. The integer flags need not be increasing
+    or contiguous and can repeat across different groups (i.e. they are
+    group-local).
 
     :arg element_flags: a :class:`numpy.ndarray` with
         :attr:`~meshmode.mesh.Mesh.nelements` entries
@@ -582,35 +593,41 @@ def split_mesh_groups(mesh, element_flags):
         are to be split.
 
     :returns: a :class:`~meshmode.mesh.Mesh` where each group has been split
-        according to flags in *element_flags*.
+        according to flags in *element_flags*. If *return_subgroup_mapping*
+        is *True*, it also returns a mapping of
+        ``(group_index, subgroup) -> new_group_index``
+
     """
     assert element_flags.shape == (mesh.nelements,)
 
     new_groups = []
-    for grp in mesh.groups:
-        mask = element_flags[
+    subgroup_to_group_map = {}
+
+    for igrp, grp in enumerate(mesh.groups):
+        grp_flags = element_flags[
                 grp.element_nr_base:grp.element_nr_base + grp.nelements]
+        unique_grp_flags = np.unique(grp_flags)
 
-        vertex_indices = grp.vertex_indices[mask, :].copy()
-        if vertex_indices.size > 0:
-            new_groups.append(grp.copy(
-                vertex_indices=vertex_indices,
-                nodes=grp.nodes[:, mask, :].copy()
-                ))
+        for iflag, flag in enumerate(unique_grp_flags):
+            subgroup_to_group_map[(igrp, flag)] = len(new_groups)
 
-        vertex_indices = grp.vertex_indices[~mask, :].copy()
-        if vertex_indices.size > 0:
+            mask = grp_flags == flag
             new_groups.append(grp.copy(
-                vertex_indices=vertex_indices,
-                nodes=grp.nodes[:, ~mask, :].copy()
+                vertex_indices=grp.vertex_indices[mask, :].copy(),
+                nodes=grp.nodes[:, mask, :].copy()
                 ))
 
     from meshmode.mesh import Mesh
-    return Mesh(
+    mesh = Mesh(
             vertices=mesh.vertices,
             groups=new_groups,
             is_conforming=mesh.is_conforming)
 
+    if return_subgroup_mapping:
+        return mesh, subgroup_to_group_map
+    else:
+        return mesh
+
 # }}}