Change tag/face mapping argument of _compute_facial_adjacency_from_vertices

meshmode/mesh/__init__.py

@@ -22,7 +22,7 @@
 
 from abc import ABC, abstractmethod
 from dataclasses import dataclass, replace, field
-from typing import Any, ClassVar, Hashable, Optional, Tuple, Type, Sequence
+from typing import Any, ClassVar, Hashable, Optional, Tuple, Type, Sequence, Mapping
 
 import numpy as np
 import numpy.linalg as la
@@ -1345,6 +1345,49 @@ def _concatenate_face_ids(face_ids_list):
         faces=np.concatenate([ids.faces for ids in face_ids_list]))
 
 
+def _find_matching_index_pairs_merged(indices):
+    """
+    Return an array containing pairs of indices into *indices* representing entries
+    that are the same.
+
+    Given an array *indices* of shape ``(N, nindices)`` containing integer-valued
+    N-tuples, returns an array of shape ``(2, nmatches)`` containing pairs of
+    indices into the second dimension of *indices* representing tuples that are
+    equal.
+
+    Returned matches are ordered such that the second element of the pair occurs
+    after the first in *indices*.
+    """
+    order = np.lexsort(indices)
+    diffs = np.diff(indices[:, order], axis=1)
+    match_indices, = (~np.any(diffs, axis=0)).nonzero()
+    # lexsort is stable, so the second entry in each match comes after the first
+    # in indices
+    return np.stack((order[match_indices], order[match_indices+1]))
+
+
+def _find_matching_index_pairs(left_indices, right_indices):
+    """
+    Return an array containing pairs of indices into *left_indices* and
+    *right_indices* representing entries that are the same.
+
+    Given an array *left_indices* of shape ``(N, nindices_left)`` and an array
+    *right_indices* of shape ``(N, nindices_right)`` containing integer-valued
+    N-tuples, returns an array of shape ``(2, nmatches)`` containing pairs of
+    indices into the second dimension of *left_indices* (first pair element) and
+    *right_indices* (second pair element) representing tuples that are equal.
+    """
+    all_indices = np.concatenate((left_indices, right_indices), axis=1)
+    index_pairs = _find_matching_index_pairs_merged(all_indices)
+
+    # Indices into left_indices are the same as indices into all_indices, but need
+    # a conversion for right_indices
+    index_pairs[1, :] -= left_indices.shape[1]
+    assert index_pairs[1, :] >= 0  # Sanity check
+
+    return index_pairs
+
+
 def _match_faces_by_vertices(groups, face_ids, vertex_index_map_func=None):
     """
     Return matching faces in *face_ids* (expressed as pairs of indices into
@@ -1388,31 +1431,28 @@ def vertex_index_map_func(vertices):
 
     # Normalize vertex-based "face identifiers" by sorting
     face_vertex_indices_increasing = np.sort(face_vertex_indices, axis=0)
-    # Lexicographically sort the face vertex indices, then diff the result to find
-    # faces with the same vertices
-    order = np.lexsort(face_vertex_indices_increasing)
-    diffs = np.diff(face_vertex_indices_increasing[:, order], axis=1)
-    match_indices, = (~np.any(diffs, axis=0)).nonzero()
 
-    return np.stack((order[match_indices], order[match_indices+1]))
+    return _find_matching_index_pairs_merged(face_vertex_indices_increasing)
 
 
 def _compute_facial_adjacency_from_vertices(
-        groups, element_id_dtype, face_id_dtype, face_vertex_indices_to_tags=None
+        groups: Sequence[MeshElementGroup],
+        element_id_dtype,
+        face_id_dtype,
+        tag_to_group_faces: Optional[Sequence[Mapping[Any, np.ndarray]]] = None
         ) -> Sequence[Sequence[FacialAdjacencyGroup]]:
+    """
+    :arg tag_to_group_faces: for each group, a mapping from tag to
+        :class:`numpy.ndarray` of shape ``(2, nfaces)`` of tagged faces in the
+        group, where ``tag_to_group_faces[igrp][tag][0, :]`` contains the element
+        indices and ``tag_to_group_faces[igrp][tag][1, :]`` contains the reference
+        face indices.
+    """
     if not groups:
         return []
 
-    if face_vertex_indices_to_tags is not None:
-        boundary_tags = {
-            tag
-            for tags in face_vertex_indices_to_tags.values()
-            for tag in tags
-            if tags is not None}
-    else:
-        boundary_tags = set()
-
-    boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)}
+    if tag_to_group_faces is None:
+        tag_to_group_faces = [{} for grp in groups]
 
     # Match up adjacent faces according to their vertex indices
 
@@ -1490,33 +1530,26 @@ def _compute_facial_adjacency_from_vertices(
             bdry_element_faces, bdry_elements = np.where(~face_has_neighbor)
             bdry_element_faces = bdry_element_faces.astype(face_id_dtype)
             bdry_elements = bdry_elements.astype(element_id_dtype)
-            belongs_to_bdry = np.full(
-                (len(boundary_tags), len(bdry_elements)), False)
-
-            if face_vertex_indices_to_tags is not None:
-                for i in range(len(bdry_elements)):
-                    ref_fvi = grp.face_vertex_indices()[bdry_element_faces[i]]
-                    fvi = frozenset(grp.vertex_indices[bdry_elements[i], ref_fvi])
-                    tags = face_vertex_indices_to_tags.get(fvi, None)
-                    if tags is not None:
-                        for tag in tags:
-                            btag_idx = boundary_tag_to_index[tag]
-                            belongs_to_bdry[btag_idx, i] = True
-
-            for btag_idx, btag in enumerate(boundary_tags):
-                indices, = np.where(belongs_to_bdry[btag_idx, :])
-                if len(indices) > 0:
-                    elements = bdry_elements[indices]
-                    element_faces = bdry_element_faces[indices]
+
+            is_tagged = np.full(len(bdry_elements), False)
+
+            for tag, tagged_elements_and_faces in tag_to_group_faces[igrp].items():
+                tagged_and_bdry_face_index_pairs = _find_matching_index_pairs(
+                    tagged_elements_and_faces.T,
+                    np.stack((bdry_elements, bdry_element_faces)))
+                bdry_face_indices = tagged_and_bdry_face_index_pairs[1, :]
+                if len(bdry_face_indices) > 0:
+                    elements = bdry_elements[bdry_face_indices]
+                    element_faces = bdry_element_faces[bdry_face_indices]
                     grp_list.append(
                         BoundaryAdjacencyGroup(
                             igroup=igrp,
-                            boundary_tag=btag,
+                            boundary_tag=tag,
                             elements=elements,
                             element_faces=element_faces))
+                    is_tagged[bdry_face_indices] = True
 
-            is_untagged = ~np.any(belongs_to_bdry, axis=0)
-            if np.any(is_untagged):
+            if not np.all(is_tagged):
                 grp_list.append(
                     BoundaryAdjacencyGroup(
                         igroup=igrp,

meshmode/mesh/generation.py

@@ -1244,63 +1244,62 @@ def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64,
 
     # {{{ compute facial adjacency for mesh if there is tag information
 
-    facial_adjacency_groups = None
-    face_vertex_indices_to_tags = {}
-    boundary_tags = list(boundary_tag_to_face.keys())
-    nbnd_tags = len(boundary_tags)
-
-    if nbnd_tags > 0:
-        vert_index_to_tuple = {
-                vertex_indices[itup]: itup
-                for itup in np.ndindex(shape)}
-
-        for ielem in range(0, grp.nelements):
-            for ref_fvi in grp.face_vertex_indices():
-                fvi = grp.vertex_indices[ielem, ref_fvi]
-                try:
-                    fvi_tuples = [vert_index_to_tuple[i] for i in fvi]
-                except KeyError:
-                    # Happens for interior faces of "X" meshes because
-                    # midpoints aren't in vert_index_to_tuple. We don't
-                    # care about them.
-                    continue
-
-                for tag in boundary_tags:
-                    # Need to map the correct face vertices to the boundary tags
-                    for face in boundary_tag_to_face[tag]:
-                        if len(face) != 2:
-                            raise ValueError(
-                                "face identifier '%s' does not "
-                                "consist of exactly two characters" % face)
-
-                        side, axis = face
-                        try:
-                            axis = axes.index(axis)
-                        except ValueError as exc:
-                            raise ValueError(
-                                "unrecognized axis in face identifier "
-                                f"'{face}'") from exc
-                        if axis >= dim:
-                            raise ValueError("axis in face identifier '%s' "
-                                             "does not exist in %dD" % (face, dim))
-
-                        if side == "-":
-                            vert_crit = 0
-                        elif side == "+":
-                            vert_crit = shape[axis] - 1
-                        else:
-                            raise ValueError("first character of face identifier"
-                                             " '%s' is not '+' or '-'" % face)
-
-                        if all(fvi_tuple[axis] == vert_crit
-                               for fvi_tuple in fvi_tuples):
-                            key = frozenset(fvi)
-                            face_vertex_indices_to_tags.setdefault(key,
-                                                                   []).append(tag)
+    tag_to_faces = {}
+
+    for tag in boundary_tag_to_face:
+        # Need to map the correct face vertices to the boundary tags
+        for face in boundary_tag_to_face[tag]:
+            if len(face) != 2:
+                raise ValueError("face identifier '%s' does not "
+                        "consist of exactly two characters" % face)
+
+            side, axis = face
+            try:
+                axis = axes.index(axis)
+            except ValueError as exc:
+                raise ValueError(
+                        f"unrecognized axis in face identifier '{face}'") from exc
+            if axis >= dim:
+                raise ValueError("axis in face identifier '%s' does not exist in %dD"
+                        % (face, dim))
+
+            if side == "-":
+                vert_crit = 0
+            elif side == "+":
+                vert_crit = shape[axis] - 1
+            else:
+                raise ValueError("first character of face identifier '%s' is not"
+                        "'+' or '-'" % face)
+
+            for fid, ref_fvi in enumerate(grp.face_vertex_indices()):
+                fvis = grp.vertex_indices[:, ref_fvi]
+                # Only consider faces whose vertices were all vertices of the
+                # original box (no midpoints, etc.)
+                candidate_elements, = np.where(
+                    np.all(fvis < nvertices, axis=1))
+                candidate_fvi_tuples = np.stack(
+                    np.unravel_index(
+                        fvis[candidate_elements, :], shape),
+                    axis=-1)
+                boundary_elements = candidate_elements[
+                    np.all(
+                        candidate_fvi_tuples[:, :, axis] == vert_crit, axis=1)]
+                boundary_faces = np.stack(
+                    (
+                        boundary_elements,
+                        np.full(boundary_elements.shape, fid)),
+                    axis=-1)
+                if tag in tag_to_faces:
+                    tag_to_faces[tag] = np.concatenate((
+                        tag_to_faces[tag],
+                        boundary_faces))
+                else:
+                    tag_to_faces[tag] = boundary_faces
 
+    if tag_to_faces:
         from meshmode.mesh import _compute_facial_adjacency_from_vertices
         facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
-                [grp], np.int32, np.int8, face_vertex_indices_to_tags)
+                [grp], np.int32, np.int8, [tag_to_faces])
     else:
         facial_adjacency_groups = None