From bd207ab9f460c9498448af2ed8604be254f6538a Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Tue, 22 Jun 2021 09:42:55 -0500
Subject: [PATCH 01/20] add mesh boundary gluing

---
 meshmode/mesh/processing.py | 244 ++++++++++++++++++++++++++++++++++++
 1 file changed, 244 insertions(+)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 016bdcdb6..0c166863c 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -45,6 +45,7 @@
 .. autofunction:: find_bounding_box
 .. autofunction:: merge_disjoint_meshes
 .. autofunction:: split_mesh_groups
+.. autofunction:: glue_mesh_boundaries
 
 .. autofunction:: map_mesh
 .. autofunction:: affine_map
@@ -891,6 +892,249 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
 # }}}
 
 
+# {{{ vertex matching
+
+# FIXME: This tree-based approach is probably slow; see if there's a way to do
+# something like this using numpy constructs
+def _match_vertices(
+        mesh, src_vertex_indices, tgt_vertex_indices, aff_map=None, tol=1e-12):
+    from meshmode.mesh.tools import AffineMap
+    if aff_map is None:
+        aff_map = AffineMap()
+
+    src_vertices = mesh.vertices[:, src_vertex_indices]
+    tgt_vertices = mesh.vertices[:, tgt_vertex_indices]
+
+    tgt_vertex_bboxes = np.stack((
+        tgt_vertices - tol,
+        tgt_vertices + tol))
+
+    from pytools.spatial_btree import SpatialBinaryTreeBucket
+    tree = SpatialBinaryTreeBucket(
+        np.min(tgt_vertex_bboxes[0], axis=1),
+        np.max(tgt_vertex_bboxes[1], axis=1))
+    for ivertex in range(len(tgt_vertex_indices)):
+        tree.insert(ivertex, tgt_vertex_bboxes[:, :, ivertex])
+
+    mapped_src_vertices = aff_map(src_vertices)
+
+    matched_tgt_vertices = np.full(len(src_vertex_indices), -1)
+    for ivertex in range(len(src_vertex_indices)):
+        mapped_src_vertex = mapped_src_vertices[:, ivertex]
+        matches = np.array(list(tree.generate_matches(mapped_src_vertex)))
+        match_bboxes = tgt_vertex_bboxes[:, :, matches]
+        in_bbox = np.all(
+            (mapped_src_vertex[:, np.newaxis] >= match_bboxes[0, :, :])
+            & (mapped_src_vertex[:, np.newaxis] <= match_bboxes[1, :, :]),
+            axis=0)
+        candidate_indices = matches[in_bbox]
+        if len(candidate_indices) == 0:
+            continue
+        displacement = (
+            mapped_src_vertex.reshape(-1, 1)
+            - tgt_vertices[:, candidate_indices])
+        distance_sq = np.sum(displacement**2, axis=0)
+        matched_tgt_vertices[ivertex] = (
+            tgt_vertex_indices[candidate_indices[np.argmin(distance_sq)]])
+
+    return matched_tgt_vertices
+
+# }}}
+
+
+# {{{ boundary face matching
+
+def _get_boundary_face_ids(mesh, btag):
+    from meshmode.mesh import _FaceIDs
+    face_ids_per_boundary_group = []
+    for igrp, fagrp_list in enumerate(mesh.facial_adjacency_groups):
+        matching_bdry_grps = [
+            fagrp for fagrp in fagrp_list
+            if isinstance(fagrp, BoundaryAdjacencyGroup)
+            and fagrp.boundary_tag == btag]
+        for bdry_grp in matching_bdry_grps:
+            face_ids = _FaceIDs(
+                groups=np.full(len(bdry_grp.elements), igrp),
+                elements=bdry_grp.elements,
+                faces=bdry_grp.element_faces)
+            face_ids_per_boundary_group.append(face_ids)
+
+    from meshmode.mesh import _concatenate_face_ids
+    return _concatenate_face_ids(face_ids_per_boundary_group)
+
+
+def _get_face_vertex_indices(mesh, face_ids):
+    max_face_vertices = max(
+        len(ref_fvi)
+        for grp in mesh.groups
+        for ref_fvi in grp.face_vertex_indices())
+
+    face_vertex_indices_per_group = []
+    for igrp, grp in enumerate(mesh.groups):
+        belongs_to_group = face_ids.groups == igrp
+        faces = face_ids.faces[belongs_to_group]
+        elements = face_ids.elements[belongs_to_group]
+        face_vertex_indices = np.full(
+            (len(faces), max_face_vertices), -1,
+            dtype=mesh.vertex_id_dtype)
+        for fid, ref_fvi in enumerate(grp.face_vertex_indices()):
+            is_face = faces == fid
+            face_vertex_indices[is_face, :len(ref_fvi)] = (
+                grp.vertex_indices[elements[is_face], :][:, ref_fvi])
+        face_vertex_indices_per_group.append(face_vertex_indices)
+
+    return np.stack(face_vertex_indices_per_group)
+
+
+def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
+    bdry_m_face_ids = _get_boundary_face_ids(mesh, btag_m)
+    bdry_n_face_ids = _get_boundary_face_ids(mesh, btag_n)
+
+    from pytools import single_valued
+    nfaces = single_valued((
+        len(bdry_m_face_ids.groups),
+        len(bdry_m_face_ids.elements),
+        len(bdry_m_face_ids.faces),
+        len(bdry_n_face_ids.groups),
+        len(bdry_n_face_ids.elements),
+        len(bdry_n_face_ids.faces)))
+
+    bdry_m_face_vertex_indices = _get_face_vertex_indices(mesh, bdry_m_face_ids)
+    bdry_n_face_vertex_indices = _get_face_vertex_indices(mesh, bdry_n_face_ids)
+
+    bdry_m_vertex_indices = np.unique(bdry_m_face_vertex_indices)
+    bdry_m_vertex_indices = bdry_m_vertex_indices[bdry_m_vertex_indices >= 0]
+    bdry_n_vertex_indices = np.unique(bdry_n_face_vertex_indices)
+    bdry_n_vertex_indices = bdry_n_vertex_indices[bdry_n_vertex_indices >= 0]
+
+    matched_bdry_n_vertex_indices = _match_vertices(
+        mesh, bdry_m_vertex_indices, bdry_n_vertex_indices,
+        aff_map=aff_map, tol=tol)
+
+    unmatched_bdry_m_vertex_indices = bdry_m_vertex_indices[
+        np.where(matched_bdry_n_vertex_indices < 0)[0]]
+    if len(unmatched_bdry_m_vertex_indices) > 0:
+        # TODO: What's the best way to log these?
+        # for vertex_index in unmatched_bdry_m_vertex_indices:
+        #     ...
+        raise RuntimeError("unable to match vertices between boundaries"
+            f" {btag_m} and {btag_n}")
+
+    from meshmode.mesh import _concatenate_face_ids
+    face_ids = _concatenate_face_ids([bdry_m_face_ids, bdry_n_face_ids])
+
+    max_vertex_index = max([np.max(grp.vertex_indices) for grp in mesh.groups])
+    vertex_index_map, = np.indices((max_vertex_index+1,),
+        dtype=mesh.element_id_dtype)
+    vertex_index_map[bdry_m_vertex_indices] = matched_bdry_n_vertex_indices
+
+    from meshmode.mesh import _match_faces_by_vertices
+    face_index_pairs = _match_faces_by_vertices(mesh.groups, face_ids,
+        vertex_index_map_func=lambda vs: vertex_index_map[vs])
+
+    assert face_index_pairs.shape[1] == nfaces
+
+    order = np.argsort(face_index_pairs[0, :])
+
+    from meshmode.mesh import _FaceIDs
+    return (
+        _FaceIDs(
+            groups=face_ids.groups[face_index_pairs[0, order]],
+            elements=face_ids.elements[face_index_pairs[0, order]],
+            faces=face_ids.faces[face_index_pairs[0, order]]),
+        _FaceIDs(
+            groups=face_ids.groups[face_index_pairs[1, order]],
+            elements=face_ids.elements[face_index_pairs[1, order]],
+            faces=face_ids.faces[face_index_pairs[1, order]]))
+
+# }}}
+
+
+# {{{ boundary gluing
+
+def _complete_glued_boundary_mappings(partial_glued_boundary_mappings):
+    partial_btag_pairs = set(
+        (btag_m, btag_n)
+        for btag_m, btag_n, _, _ in partial_glued_boundary_mappings)
+
+    glued_boundary_mappings = []
+    for btag_m, btag_n, aff_map, tol in partial_glued_boundary_mappings:
+        glued_boundary_mappings.append((btag_m, btag_n, aff_map, tol))
+        if (btag_n, btag_m) not in partial_btag_pairs:
+            inv_aff_map = aff_map.inverted()
+            glued_boundary_mappings.append((btag_n, btag_m, inv_aff_map, tol))
+
+    return glued_boundary_mappings
+
+
+def glue_mesh_boundaries(mesh, glued_boundary_mappings):
+    """
+    Create a new mesh from *mesh* in which one or more pairs of boundaries are
+    "glued" together such that the boundary surfaces become part of the interior
+    of the mesh. This can be used to construct, e.g., periodic boundaries.
+
+    Corresponding boundaries' vertices must map into each other via an affine
+    transformation (though the vertex ordering need not be the same).
+
+    :arg glued_boundary_mappings: a :class:`list` of tuples
+        *(btag_m, btag_n, aff_map, tol)* which each specify a mapping between two
+        boundaries in *mesh* that should be glued together. *aff_map* is a
+        :class:`~meshmode.mesh.tools.AffineMap` that represents the affine mapping
+        from the vertices of boundary *btag_m* into the vertices of boundary
+        *btag_n*. *tol* is the tolerance allowed between the vertex coordinates of
+        *btag_n* and the transformed vertex coordinates of *btag_m* when attempting
+        to match the two.
+    """
+    glued_boundary_mappings = _complete_glued_boundary_mappings(
+        glued_boundary_mappings)
+
+    glued_btags = (
+        set(btag_m for btag_m, _, _, _ in glued_boundary_mappings)
+        | set(btag_n for _, btag_n, _, _ in glued_boundary_mappings))
+
+    face_id_pairs_for_mapping = []
+    for btag_m, btag_n, aff_map, tol in glued_boundary_mappings:
+        face_id_pairs_for_mapping.append(_match_boundary_faces(mesh,
+            btag_m, btag_n, aff_map, tol))
+
+    from meshmode.mesh import InteriorAdjacencyGroup, BoundaryAdjacencyGroup
+
+    facial_adjacency_groups = []
+
+    for igrp, old_fagrp_list in enumerate(mesh.facial_adjacency_groups):
+        fagrp_list = [
+            fagrp for fagrp in old_fagrp_list
+            if not isinstance(fagrp, BoundaryAdjacencyGroup)
+            or fagrp.boundary_tag not in glued_btags]
+
+        for imap, (_, _, aff_map, _) in enumerate(glued_boundary_mappings):
+            mapping_face_id_pairs = face_id_pairs_for_mapping[imap]
+            belongs_to_group = mapping_face_id_pairs[0].groups == igrp
+            for ineighbor_grp in range(len(mesh.groups)):
+                indices, = np.where(
+                    belongs_to_group
+                    & (mapping_face_id_pairs[1].groups == ineighbor_grp))
+                if len(indices) > 0:
+                    elements = mapping_face_id_pairs[0].elements[indices]
+                    element_faces = mapping_face_id_pairs[0].faces[indices]
+                    neighbors = mapping_face_id_pairs[1].elements[indices]
+                    neighbor_faces = mapping_face_id_pairs[1].faces[indices]
+                    fagrp_list.append(InteriorAdjacencyGroup(
+                        igroup=igrp,
+                        ineighbor_group=ineighbor_grp,
+                        elements=elements,
+                        element_faces=element_faces,
+                        neighbors=neighbors,
+                        neighbor_faces=neighbor_faces,
+                        aff_map=aff_map))
+
+        facial_adjacency_groups.append(fagrp_list)
+
+    return mesh.copy(facial_adjacency_groups=facial_adjacency_groups)
+
+# }}}
+
+
 # {{{ map
 
 def map_mesh(mesh, f):  # noqa

From 2af9d3973e4f0fdcd1433777f653f1388888761b Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Tue, 22 Jun 2021 09:36:06 -0500
Subject: [PATCH 02/20] add generate_annular_cylinder_slice_mesh

---
 meshmode/mesh/generation.py | 53 +++++++++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index c978f181e..8d0c5f805 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -70,6 +70,7 @@
 .. autofunction:: generate_box_mesh
 .. autofunction:: generate_regular_rect_mesh
 .. autofunction:: generate_warped_rect_mesh
+.. autofunction:: generate_annular_cylinder_slice_mesh
 
 Tools for Iterative Refinement
 ------------------------------
@@ -1267,6 +1268,58 @@ def m(x):
 # }}}
 
 
+# {{{ generate_annular_cylinder_slice_mesh
+
+def generate_annular_cylinder_slice_mesh(
+        n, center, inner_radius, outer_radius, periodic=False):
+    r"""
+    Generate a slice of a 3D annular cylinder for
+    :math:`\theta \in [-\frac{\pi}{4}, \frac{\pi}{4}]`. Optionally periodic in
+    $\theta$.
+    """
+    unit_mesh = generate_regular_rect_mesh(
+        a=(0,)*3,
+        b=(1,)*3,
+        nelements_per_axis=(n,)*3,
+        boundary_tag_to_face={
+            "-r": ["-x"],
+            "+r": ["+x"],
+            "-theta": ["-y"],
+            "+theta": ["+y"],
+            "-z": ["-z"],
+            "+z": ["+z"],
+            })
+
+    def transform(x):
+        r = inner_radius*(1 - x[0]) + outer_radius*x[0]
+        theta = -np.pi/4*(1 - x[1]) + np.pi/4*x[1]
+        z = -0.5*(1 - x[2]) + 0.5*x[2]
+        return (
+            center[0] + r*np.cos(theta),
+            center[1] + r*np.sin(theta),
+            center[2] + z)
+
+    from meshmode.mesh.processing import map_mesh
+    mesh = map_mesh(unit_mesh, lambda x: np.stack(transform(x)))
+
+    if periodic:
+        from meshmode.mesh.processing import _get_rotation_matrix_from_angle_and_axis
+        matrix = _get_rotation_matrix_from_angle_and_axis(
+            np.pi/2, np.array([0, 0, 1]))
+        from meshmode.mesh.tools import AffineMap
+        aff_map = AffineMap(matrix, center - matrix @ center)
+
+        from meshmode.mesh.processing import glue_mesh_boundaries
+        periodic_mesh = glue_mesh_boundaries(mesh,
+            glued_boundary_mappings=[("-theta", "+theta", aff_map, 1e-12)])
+
+        return periodic_mesh
+    else:
+        return mesh
+
+# }}}
+
+
 # {{{ warp_and_refine_until_resolved
 
 @log_process(logger)

From 8d645a38137c8a9a50a0c32fe9395c2aae831e0a Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 9 Aug 2021 11:02:46 -0500
Subject: [PATCH 03/20] add tests for facial adjacency transforms and mesh
 boundary gluing

---
 test/test_mesh.py     | 142 +++++++++++++++++++++++++++++++++++++++++-
 test/test_meshmode.py |  14 ++++-
 2 files changed, 152 insertions(+), 4 deletions(-)

diff --git a/test/test_mesh.py b/test/test_mesh.py
index 162eb1e7c..921c1f25e 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -35,18 +35,37 @@
     Mesh,
     SimplexElementGroup,
     TensorProductElementGroup,
+    InteriorAdjacencyGroup,
     BoundaryAdjacencyGroup)
 from meshmode.discretization.poly_element import (
         default_simplex_group_factory,
         LegendreGaussLobattoTensorProductGroupFactory,
         )
 import meshmode.mesh.generation as mgen
+from meshmode.mesh.tools import AffineMap
 
 
 import logging
 logger = logging.getLogger(__name__)
 
 
+def _get_rotation(amount, axis, center=None):
+    """
+    Return a matrix (if *center* is ``None``) or
+    :class:`~meshmode.mesh.tools.AffineMap` (if *center* is not ``None``)
+    corresponding to a rotation by *amount* (in radians) through a vector *axis*
+    centered at *center*. *center* defaults to the origin if not specified.
+    """
+    from meshmode.mesh.processing import _get_rotation_matrix_from_angle_and_axis
+    matrix = _get_rotation_matrix_from_angle_and_axis(amount, axis)
+    if center is None:
+        return matrix
+    else:
+        # x0 + matrix @ (x - x0) = matrix @ x + (I - matrix) @ x0
+        offset = (np.eye(3) - matrix) @ center
+        return AffineMap(matrix, offset)
+
+
 # {{{ test_nonequal_rect_mesh_generation
 
 @pytest.mark.parametrize(("dim", "mesh_type"), [
@@ -170,7 +189,6 @@ def test_mesh_as_python():
 # {{{ test_affine_map
 
 def test_affine_map():
-    from meshmode.mesh.tools import AffineMap
     for d in range(1, 5):
         for _ in range(100):
             a = np.random.randn(d, d)+10*np.eye(d)
@@ -778,6 +796,128 @@ def test_cube_icosphere(actx_factory, order, visualize=True):
 # }}}
 
 
+# {{{ mesh boundary gluing
+
+def test_glued_mesh():
+    n = 4
+    center = (1, 2, 3)
+
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(n, center, 0.5, 1)
+
+    map_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]), center)
+    map_upper_to_lower = _get_rotation(-np.pi/2, np.array([0, 0, 1]), center)
+
+    from meshmode.mesh.processing import glue_mesh_boundaries
+    mesh = glue_mesh_boundaries(orig_mesh,
+        glued_boundary_mappings=[
+            ("-theta", "+theta", map_lower_to_upper, 1e-12)
+        ])
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+    assert len(int_grps) == 3
+
+    bdry_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, BoundaryAdjacencyGroup)]
+    assert len(bdry_grps) == (
+        4  # +/-r and +/-z
+        + 3  # BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL
+        )
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    from pytools import single_valued
+
+    n_lower_faces = single_valued((
+        len(lower_grp.elements),
+        len(lower_grp.element_faces),
+        len(lower_grp.neighbors),
+        len(lower_grp.neighbor_faces)))
+    assert n_lower_faces == 2*n**2
+    assert lower_grp.aff_map == map_lower_to_upper
+
+    n_upper_faces = single_valued((
+        len(upper_grp.elements),
+        len(upper_grp.element_faces),
+        len(upper_grp.neighbors),
+        len(upper_grp.neighbor_faces)))
+    assert n_upper_faces == 2*n**2
+    assert upper_grp.aff_map == map_upper_to_lower
+
+    lower_face_indices = np.full(
+        (orig_mesh.groups[0].nfaces, orig_mesh.groups[0].nelements), -1)
+    upper_face_indices = np.full(
+        (orig_mesh.groups[0].nfaces, orig_mesh.groups[0].nelements), -1)
+
+    lower_face_indices[lower_grp.element_faces, lower_grp.elements] = (
+        np.indices((n_lower_faces,)))
+    upper_face_indices[upper_grp.element_faces, upper_grp.elements] = (
+        np.indices((n_upper_faces,)))
+
+    indices = upper_face_indices[lower_grp.neighbor_faces, lower_grp.neighbors]
+    assert np.all(indices >= 0)
+    assert np.all(upper_grp.neighbors[indices] == lower_grp.elements)
+    assert np.all(upper_grp.neighbor_faces[indices] == lower_grp.element_faces)
+
+    indices = lower_face_indices[upper_grp.neighbor_faces, upper_grp.neighbors]
+    assert np.all(indices >= 0)
+    assert np.all(lower_grp.neighbors[indices] == upper_grp.elements)
+    assert np.all(lower_grp.neighbor_faces[indices] == upper_grp.element_faces)
+
+
+def test_glued_mesh_matrix_only():
+    n = 4
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(n, (0, 0, 0), 0.5, 1)
+
+    matrix_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]))
+    matrix_upper_to_lower = _get_rotation(-np.pi/2, np.array([0, 0, 1]))
+
+    from meshmode.mesh.processing import glue_mesh_boundaries
+    mesh = glue_mesh_boundaries(orig_mesh,
+        glued_boundary_mappings=[
+            ("-theta", "+theta", AffineMap(matrix=matrix_lower_to_upper), 1e-12)
+        ])
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    assert lower_grp.aff_map == AffineMap(matrix=matrix_lower_to_upper)
+    assert upper_grp.aff_map == AffineMap(matrix=matrix_upper_to_lower)
+
+
+def test_glued_mesh_offset_only():
+    n = 4
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(n, (0, 0, 0), 0.5, 1)
+
+    offset_lower_to_upper = np.array([0, 0, 1])
+    offset_upper_to_lower = np.array([0, 0, -1])
+
+    from meshmode.mesh.processing import glue_mesh_boundaries
+    mesh = glue_mesh_boundaries(orig_mesh,
+        glued_boundary_mappings=[
+            ("-z", "+z", AffineMap(offset=offset_lower_to_upper), 1e-12)
+        ])
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    assert lower_grp.aff_map == AffineMap(offset=offset_lower_to_upper)
+    assert upper_grp.aff_map == AffineMap(offset=offset_upper_to_lower)
+
+# }}}
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1:
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index ceb7b9cef..3d28b6468 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -334,13 +334,14 @@ def f(x):
     ("segment", 1, [8, 16, 32]),
     ("blob", 2, [1e-1, 8e-2, 5e-2]),
     ("warp", 2, [3, 5, 7]),
-    ("warp", 3, [5, 7])
+    ("warp", 3, [5, 7]),
+    ("periodic", 3, [5, 7])
     ])
 def test_opposite_face_interpolation(actx_factory, group_factory,
         mesh_name, dim, mesh_pars):
     if (group_factory is LegendreGaussLobattoTensorProductGroupFactory
-            and mesh_name in ["segment", "blob"]):
-        pytest.skip("tensor products not implemented on blobs")
+            and mesh_name in ["segment", "blob", "periodic"]):
+        pytest.skip(f"tensor products not implemented on {mesh_name}")
 
     logging.basicConfig(level=logging.INFO)
     actx = actx_factory()
@@ -395,6 +396,13 @@ def f(x):
             mesh = mgen.generate_warped_rect_mesh(dim, order=order,
                     nelements_side=mesh_par, group_cls=group_cls)
 
+            h = 1/mesh_par
+        elif mesh_name == "periodic":
+            assert dim == 3
+
+            mesh = mgen.generate_annular_cylinder_slice_mesh(
+                mesh_par, (1, 2, 3), 0.5, 1, periodic=True)
+
             h = 1/mesh_par
         else:
             raise ValueError("mesh_name not recognized")

From d51a69193f98c1ddcdf9e296791df1528c2f0ce9 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Thu, 17 Jun 2021 09:18:01 -0500
Subject: [PATCH 04/20] check for facial adjacency transforms in map_mesh

---
 meshmode/mesh/processing.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 0c166863c..68185cff1 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1141,6 +1141,18 @@ def map_mesh(mesh, f):  # noqa
     """Apply the map *f* to the mesh. *f* needs to accept and return arrays of
     shape ``(ambient_dim, npoints)``."""
 
+    if mesh._facial_adjacency_groups is not None:
+        has_adj_maps = any([
+            hasattr(fagrp, "aff_map")
+            and (fagrp.aff_map.matrix is not None
+                or fagrp.aff_map.offset is not None)
+            for fagrp_list in mesh.facial_adjacency_groups
+            for fagrp in fagrp_list])
+        if has_adj_maps:
+            raise ValueError("cannot apply a general map to a mesh that has "
+                "affine mappings in its facial adjacency. If the map is affine, "
+                "use affine_map instead")
+
     vertices = f(mesh.vertices)
     if not vertices.flags.c_contiguous:
         vertices = np.copy(vertices, order="C")

From 607755a15adb6673af4de3dd3c2c9d5a9ac89c34 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 9 Aug 2021 11:03:28 -0500
Subject: [PATCH 05/20] add facial adjacency transform handling in affine_map

---
 meshmode/mesh/processing.py | 67 +++++++++++++++++++++++++-
 test/test_mesh.py           | 95 +++++++++++++++++++++++++++++++++++++
 2 files changed, 161 insertions(+), 1 deletion(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 68185cff1..cc584209a 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1201,7 +1201,72 @@ def affine_map(mesh,
         raise ValueError(f"b has shape '{b.shape}' for a {mesh.ambient_dim}d mesh")
 
     from meshmode.mesh.tools import AffineMap
-    return map_mesh(mesh, AffineMap(A, b))
+    f = AffineMap(A, b)
+
+    vertices = f(mesh.vertices)
+    if not vertices.flags.c_contiguous:
+        vertices = np.copy(vertices, order="C")
+
+    # {{{ assemble new groups list
+
+    new_groups = []
+
+    for group in mesh.groups:
+        mapped_nodes = f(group.nodes.reshape(mesh.ambient_dim, -1))
+        if not mapped_nodes.flags.c_contiguous:
+            mapped_nodes = np.copy(mapped_nodes, order="C")
+
+        new_groups.append(group.copy(
+            nodes=mapped_nodes.reshape(*group.nodes.shape)))
+
+    # }}}
+
+    # {{{ assemble new facial adjacency groups
+
+    if mesh._facial_adjacency_groups is not None:
+        # For a facial adjacency transform T(x) = Gx + h in the original mesh,
+        # its corresponding transform in the new mesh will be (T')(x) = G'x + h',
+        # where:
+        # G' = G
+        # h' = Ah + (I - G)b
+        def compute_new_map(old_map):
+            if old_map.matrix is not None:
+                matrix = old_map.matrix.copy()
+            else:
+                matrix = None
+            if old_map.offset is not None:
+                if A is not None:
+                    offset = A @ old_map.offset
+                else:
+                    offset = old_map.offset.copy()
+                if matrix is not None and b is not None:
+                    offset += b - matrix @ b
+            else:
+                offset = None
+            return AffineMap(matrix, offset)
+
+        facial_adjacency_groups = []
+        for old_fagrp_list in mesh.facial_adjacency_groups:
+            fagrp_list = []
+            for old_fagrp in old_fagrp_list:
+                if hasattr(old_fagrp, "aff_map"):
+                    aff_map = compute_new_map(old_fagrp.aff_map)
+                    fagrp_list.append(
+                        old_fagrp.copy(
+                            aff_map=aff_map))
+                else:
+                    fagrp_list.append(old_fagrp.copy())
+            facial_adjacency_groups.append(fagrp_list)
+
+    else:
+        facial_adjacency_groups = None
+
+    # }}}
+
+    return mesh.copy(
+            vertices=vertices, groups=new_groups,
+            facial_adjacency_groups=facial_adjacency_groups,
+            is_conforming=mesh.is_conforming)
 
 
 def _get_rotation_matrix_from_angle_and_axis(theta, axis):
diff --git a/test/test_mesh.py b/test/test_mesh.py
index 921c1f25e..672f7801b 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -223,6 +223,101 @@ def test_partial_affine_map(dim=2):
     assert la.norm(orig_mesh.vertices - mesh.vertices / np.pi) < 1.0e-14
 
 
+def test_affine_map_with_facial_adjacency_maps(visualize=False):
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(
+        4, (1, 2, 0), 0.5, 1, periodic=True)
+
+    if visualize:
+        from meshmode.mesh.visualization import write_vertex_vtk_file
+        write_vertex_vtk_file(orig_mesh, "affine_map_facial_adj_original.vtu")
+
+    from meshmode.mesh.processing import affine_map
+
+    tol = 1e-12
+
+    def almost_equal(map1, map2):
+        def component_almost_equal(array1, array2):
+            if isinstance(array1, np.ndarray) and isinstance(array2, np.ndarray):
+                return la.norm(array1 - array2) < tol
+            else:
+                return array1 == array2
+
+        return (
+            component_almost_equal(map1.matrix, map2.matrix)
+            and component_almost_equal(map1.offset, map2.offset))
+
+    # Matrix only
+    mesh = affine_map(orig_mesh, A=_get_rotation(np.pi/2, axis=np.array([0, 0, 1])))
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_matrix.vtu")
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+    assert len(int_grps) == 3
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    assert almost_equal(
+        lower_grp.aff_map,
+        _get_rotation(
+            np.pi/2, axis=np.array([0, 0, 1]), center=np.array([-2, 1, 0])))
+    assert almost_equal(
+        upper_grp.aff_map,
+        _get_rotation(
+            -np.pi/2, axis=np.array([0, 0, 1]), center=np.array([-2, 1, 0])))
+
+    # Offset only
+    mesh = affine_map(orig_mesh, b=np.array([0, -2, 0]))
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_offset.vtu")
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+    assert len(int_grps) == 3
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    assert almost_equal(
+        lower_grp.aff_map,
+        _get_rotation(
+            np.pi/2, axis=np.array([0, 0, 1]), center=np.array([1, 0, 0])))
+    assert almost_equal(
+        upper_grp.aff_map,
+        _get_rotation(
+            -np.pi/2, axis=np.array([0, 0, 1]), center=np.array([1, 0, 0])))
+
+    # Matrix and offset
+    aff_map = _get_rotation(
+        np.pi/2, axis=np.array([0, 0, 1]), center=np.array([1, 1, 0]))
+    mesh = affine_map(orig_mesh, A=aff_map.matrix, b=aff_map.offset)
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_matrix_and_offset.vtu")
+
+    int_grps = [
+        fagrp for fagrp in mesh.facial_adjacency_groups[0]
+        if isinstance(fagrp, InteriorAdjacencyGroup)]
+    assert len(int_grps) == 3
+
+    lower_grp = int_grps[1]
+    upper_grp = int_grps[2]
+
+    assert almost_equal(
+        lower_grp.aff_map,
+        _get_rotation(
+            np.pi/2, axis=np.array([0, 0, 1]), center=np.array([0, 1, 0])))
+    assert almost_equal(
+        upper_grp.aff_map,
+        _get_rotation(
+            -np.pi/2, axis=np.array([0, 0, 1]), center=np.array([0, 1, 0])))
+
+
 @pytest.mark.parametrize("ambient_dim", [2, 3])
 def test_mesh_rotation(ambient_dim, visualize=False):
     order = 3

From b18b5a7bb63ccbee4d4a0bd31b53dc40f02882d7 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 9 Aug 2021 10:37:51 -0500
Subject: [PATCH 06/20] add periodic argument to
 generate_box_mesh/generate_regular_rect_mesh

---
 meshmode/mesh/generation.py | 45 ++++++++++++++++++++++++++++++++++---
 test/test_meshmode.py       | 18 +++++++++++----
 2 files changed, 56 insertions(+), 7 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 8d0c5f805..4aaa4c7f7 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -882,13 +882,15 @@ def generate_urchin(
 
 @deprecate_keyword("group_factory", "group_cls")
 def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
-        group_cls=None, boundary_tag_to_face=None,
+        periodic=None, group_cls=None, boundary_tag_to_face=None,
         mesh_type=None, unit_nodes=None):
     r"""Create a semi-structured mesh.
 
     :param axis_coords: a tuple with a number of entries corresponding
         to the number of dimensions, with each entry a numpy array
         specifying the coordinates to be used along that axis.
+    :param periodic: an optional tuple of :class:`bool` indicating whether
+        the mesh is periodic along each axis.
     :param group_cls: One of :class:`meshmode.mesh.SimplexElementGroup`
         or :class:`meshmode.mesh.TensorProductElementGroup`.
     :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
@@ -947,6 +949,9 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     dim = len(axis_coords)
 
+    if periodic is None:
+        periodic = (False,)*dim
+
     shape = tuple(len(axc) for axc in axis_coords)
 
     from pytools import product
@@ -1079,12 +1084,20 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
             vertices.reshape(dim, -1), el_vertices, order,
             group_cls=group_cls, unit_nodes=unit_nodes)
 
+    axes = ["x", "y", "z", "w"]
+
+    for idim in range(dim):
+        if periodic[idim]:
+            lower_face = "-" + axes[idim]
+            upper_face = "+" + axes[idim]
+            boundary_tag_to_face["periodic_" + lower_face] = [lower_face]
+            boundary_tag_to_face["periodic_" + upper_face] = [upper_face]
+
     # {{{ 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())
-    axes = ["x", "y", "z", "w"]
 
     if boundary_tags:
         vert_index_to_tuple = {
@@ -1141,10 +1154,32 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     # }}}
 
     from meshmode.mesh import Mesh
-    return Mesh(vertices, [grp],
+    mesh = Mesh(vertices, [grp],
             facial_adjacency_groups=facial_adjacency_groups,
             is_conforming=True)
 
+    if any(periodic):
+        from meshmode.mesh.tools import AffineMap
+        glued_boundary_mappings = []
+        for idim in range(dim):
+            if periodic[idim]:
+                lower_face = "-" + axes[idim]
+                upper_face = "+" + axes[idim]
+                offset = np.zeros(dim, dtype=np.float64)
+                offset[idim] = axis_coords[idim][-1] - axis_coords[idim][0]
+                glued_boundary_mappings.append((
+                    "periodic_" + lower_face,
+                    "periodic_" + upper_face,
+                    AffineMap(offset=offset),
+                    1e-12*offset[idim]))
+
+        from meshmode.mesh.processing import glue_mesh_boundaries
+        periodic_mesh = glue_mesh_boundaries(mesh, glued_boundary_mappings)
+
+        return periodic_mesh
+    else:
+        return mesh
+
 # }}}
 
 
@@ -1153,6 +1188,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 @deprecate_keyword("group_factory", "group_cls")
 def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), *, nelements_per_axis=None,
                                npoints_per_axis=None,
+                               periodic=None,
                                order=1,
                                boundary_tag_to_face=None,
                                group_cls=None,
@@ -1167,6 +1203,8 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), *, nelements_per_axis=None,
         number of elements along each axis.
     :param npoints_per_axis: an optional tuple of integers indicating the
         number of points along each axis.
+    :param periodic: an optional tuple of :class:`bool` indicating whether
+        the mesh is periodic along each axis.
     :param order: the mesh element order.
     :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
         See :func:`generate_box_mesh`.
@@ -1206,6 +1244,7 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), *, nelements_per_axis=None,
             for a_i, b_i, npoints_i in zip(a, b, npoints_per_axis)]
 
     return generate_box_mesh(axis_coords, order=order,
+                             periodic=periodic,
                              boundary_tag_to_face=boundary_tag_to_face,
                              group_cls=group_cls,
                              mesh_type=mesh_type)
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 3d28b6468..3a80c8512 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -335,6 +335,7 @@ def f(x):
     ("blob", 2, [1e-1, 8e-2, 5e-2]),
     ("warp", 2, [3, 5, 7]),
     ("warp", 3, [5, 7]),
+    ("periodic", 2, [3, 5, 7]),
     ("periodic", 3, [5, 7])
     ])
 def test_opposite_face_interpolation(actx_factory, group_factory,
@@ -398,12 +399,21 @@ def f(x):
 
             h = 1/mesh_par
         elif mesh_name == "periodic":
-            assert dim == 3
+            assert dim == 2 or dim == 3
 
-            mesh = mgen.generate_annular_cylinder_slice_mesh(
-                mesh_par, (1, 2, 3), 0.5, 1, periodic=True)
+            if dim == 2:
+                mesh = mgen.generate_regular_rect_mesh(
+                    a=(-np.pi/2,)*dim,
+                    b=((3*np.pi)/2,)*dim,
+                    nelements_per_axis=(mesh_par,)*dim,
+                    periodic=(True, False))
 
-            h = 1/mesh_par
+                h = 1/mesh_par
+            else:
+                mesh = mgen.generate_annular_cylinder_slice_mesh(
+                    mesh_par, (1, 2, 3), 0.5, 1, periodic=True)
+
+                h = 1/mesh_par
         else:
             raise ValueError("mesh_name not recognized")
 

From d772f4f1f2274987e0e851ca87f8f305092fb5f7 Mon Sep 17 00:00:00 2001
From: Matt Smith <mjsmith6@illinois.edu>
Date: Mon, 4 Oct 2021 13:37:07 -0500
Subject: [PATCH 07/20] make arguments keyword-only in generate_box_mesh
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Andreas Klöckner <inform@tiker.net>
---
 meshmode/mesh/generation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 4aaa4c7f7..c86fe495c 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -881,7 +881,7 @@ def generate_urchin(
 # {{{ generate_box_mesh
 
 @deprecate_keyword("group_factory", "group_cls")
-def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
+def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64,
         periodic=None, group_cls=None, boundary_tag_to_face=None,
         mesh_type=None, unit_nodes=None):
     r"""Create a semi-structured mesh.

From b012ba3dfb1947c1e7c3019f69e30c2c244e342b Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 4 Oct 2021 14:06:11 -0500
Subject: [PATCH 08/20] reference glue_mesh_boundaries in generate_box_mesh and
 generate_regular_rect_mesh

---
 meshmode/mesh/generation.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index ff0456470..ab54c97c3 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -890,7 +890,8 @@ def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64,
         to the number of dimensions, with each entry a numpy array
         specifying the coordinates to be used along that axis.
     :arg periodic: an optional tuple of :class:`bool` indicating whether
-        the mesh is periodic along each axis.
+        the mesh is periodic along each axis. Acts as a shortcut for calling
+        :func:`meshmode.mesh.processing.glue_mesh_boundaries`.
     :arg group_cls: One of :class:`meshmode.mesh.SimplexElementGroup`
         or :class:`meshmode.mesh.TensorProductElementGroup`.
     :arg boundary_tag_to_face: an optional dictionary for tagging boundaries.
@@ -1204,7 +1205,8 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), *, nelements_per_axis=None,
     :arg npoints_per_axis: an optional tuple of integers indicating the
         number of points along each axis.
     :arg periodic: an optional tuple of :class:`bool` indicating whether
-        the mesh is periodic along each axis.
+        the mesh is periodic along each axis. Acts as a shortcut for calling
+        :func:`meshmode.mesh.processing.glue_mesh_boundaries`.
     :arg order: the mesh element order.
     :arg boundary_tag_to_face: an optional dictionary for tagging boundaries.
         See :func:`generate_box_mesh`.

From 8448c94072a743a4b4c5bfe451f4723e555ba8a2 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 4 Oct 2021 14:24:42 -0500
Subject: [PATCH 09/20] fix doc reference

---
 meshmode/mesh/processing.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 006367b87..79b6d7a3b 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1076,11 +1076,10 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
     :arg glued_boundary_mappings: a :class:`list` of tuples
         *(btag_m, btag_n, aff_map, tol)* which each specify a mapping between two
         boundaries in *mesh* that should be glued together. *aff_map* is a
-        :class:`~meshmode.mesh.tools.AffineMap` that represents the affine mapping
-        from the vertices of boundary *btag_m* into the vertices of boundary
-        *btag_n*. *tol* is the tolerance allowed between the vertex coordinates of
-        *btag_n* and the transformed vertex coordinates of *btag_m* when attempting
-        to match the two.
+        :class:`~meshmode.AffineMap` that represents the affine mapping from the
+        vertices of boundary *btag_m* into the vertices of boundary *btag_n*. *tol*
+        is the tolerance allowed between the vertex coordinates of *btag_n* and the
+        transformed vertex coordinates of *btag_m* when attempting to match the two.
     """
     glued_boundary_mappings = _complete_glued_boundary_mappings(
         glued_boundary_mappings)

From 8d08353f6e7c617f8f68c8bb4c1625ad723ac7a5 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Mon, 4 Oct 2021 16:51:22 -0500
Subject: [PATCH 10/20] report information about vertices that could not be
 matched

---
 meshmode/mesh/processing.py | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 79b6d7a3b..85b70ac64 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1010,12 +1010,17 @@ def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
 
     unmatched_bdry_m_vertex_indices = bdry_m_vertex_indices[
         np.where(matched_bdry_n_vertex_indices < 0)[0]]
-    if len(unmatched_bdry_m_vertex_indices) > 0:
-        # TODO: What's the best way to log these?
-        # for vertex_index in unmatched_bdry_m_vertex_indices:
-        #     ...
-        raise RuntimeError("unable to match vertices between boundaries"
-            f" {btag_m} and {btag_n}")
+    nunmatched = len(unmatched_bdry_m_vertex_indices)
+    if nunmatched > 0:
+        vertices = mesh.vertices[:, unmatched_bdry_m_vertex_indices]
+        mapped_vertices = aff_map(vertices)
+        raise RuntimeError(
+            f"unable to match vertices between boundaries {btag_m} and {btag_n}.\n"
+            + "Unmatched vertices (original -> mapped):\n"
+            + "\n".join([
+                f"{vertices[:, i]} -> {mapped_vertices[:, i]}"
+                for i in range(min(nunmatched, 10))])
+            + f"\n...\n({nunmatched-10} more omitted.)" if nunmatched > 10 else "")
 
     from meshmode.mesh import _concatenate_face_ids
     face_ids = _concatenate_face_ids([bdry_m_face_ids, bdry_n_face_ids])

From 4f799eaa245c8bd8c36078b3144de7b654ea7818 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Tue, 5 Oct 2021 09:19:54 -0500
Subject: [PATCH 11/20] disable nodal adjacency in glued meshes

---
 meshmode/mesh/processing.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 85b70ac64..8603ad03e 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1076,7 +1076,9 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
     of the mesh. This can be used to construct, e.g., periodic boundaries.
 
     Corresponding boundaries' vertices must map into each other via an affine
-    transformation (though the vertex ordering need not be the same).
+    transformation (though the vertex ordering need not be the same). Currently
+    operates only on facial adjacency; any existing nodal adjacency in *mesh* is
+    ignored/invalidated.
 
     :arg glued_boundary_mappings: a :class:`list` of tuples
         *(btag_m, btag_n, aff_map, tol)* which each specify a mapping between two
@@ -1131,7 +1133,9 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
 
         facial_adjacency_groups.append(fagrp_list)
 
-    return mesh.copy(facial_adjacency_groups=facial_adjacency_groups)
+    return mesh.copy(
+        nodal_adjacency=False,
+        facial_adjacency_groups=facial_adjacency_groups)
 
 # }}}
 

From 6484478951733ccf623b90f71e289994b9d9e617 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Tue, 5 Oct 2021 09:37:58 -0500
Subject: [PATCH 12/20] use list/set comprehensions

---
 meshmode/mesh/processing.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 8603ad03e..61d6c9e99 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1055,9 +1055,9 @@ def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
 # {{{ boundary gluing
 
 def _complete_glued_boundary_mappings(partial_glued_boundary_mappings):
-    partial_btag_pairs = set(
+    partial_btag_pairs = {
         (btag_m, btag_n)
-        for btag_m, btag_n, _, _ in partial_glued_boundary_mappings)
+        for btag_m, btag_n, _, _ in partial_glued_boundary_mappings}
 
     glued_boundary_mappings = []
     for btag_m, btag_n, aff_map, tol in partial_glued_boundary_mappings:
@@ -1091,14 +1091,14 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
     glued_boundary_mappings = _complete_glued_boundary_mappings(
         glued_boundary_mappings)
 
-    glued_btags = (
-        set(btag_m for btag_m, _, _, _ in glued_boundary_mappings)
-        | set(btag_n for _, btag_n, _, _ in glued_boundary_mappings))
+    glued_btags = {
+        btag
+        for btag_m, btag_n, _, _ in glued_boundary_mappings
+        for btag in (btag_m, btag_n)}
 
-    face_id_pairs_for_mapping = []
-    for btag_m, btag_n, aff_map, tol in glued_boundary_mappings:
-        face_id_pairs_for_mapping.append(_match_boundary_faces(mesh,
-            btag_m, btag_n, aff_map, tol))
+    face_id_pairs_for_mapping = [
+        _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol)
+        for btag_m, btag_n, aff_map, tol in glued_boundary_mappings]
 
     from meshmode.mesh import InteriorAdjacencyGroup, BoundaryAdjacencyGroup
 

From 466457541d539b9056c91cd7f82169a604abdbec Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Tue, 5 Oct 2021 14:33:11 -0500
Subject: [PATCH 13/20] use dataclass for boundary mappings

---
 meshmode/mesh/generation.py | 29 ++++++-----
 meshmode/mesh/processing.py | 96 +++++++++++++++++++++++++------------
 test/test_mesh.py           | 41 +++++++++-------
 3 files changed, 106 insertions(+), 60 deletions(-)

diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index ab54c97c3..faf7a27b9 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -1160,22 +1160,23 @@ def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64,
             is_conforming=True)
 
     if any(periodic):
-        from meshmode.mesh.tools import AffineMap
-        glued_boundary_mappings = []
+        from meshmode.mesh.processing import (
+            glue_mesh_boundaries, BoundaryPairMapping)
+
+        from meshmode import AffineMap
+        bdry_pair_mappings_and_tols = []
         for idim in range(dim):
             if periodic[idim]:
-                lower_face = "-" + axes[idim]
-                upper_face = "+" + axes[idim]
                 offset = np.zeros(dim, dtype=np.float64)
                 offset[idim] = axis_coords[idim][-1] - axis_coords[idim][0]
-                glued_boundary_mappings.append((
-                    "periodic_" + lower_face,
-                    "periodic_" + upper_face,
-                    AffineMap(offset=offset),
+                bdry_pair_mappings_and_tols.append((
+                    BoundaryPairMapping(
+                        "periodic_-" + axes[idim],
+                        "periodic_+" + axes[idim],
+                        AffineMap(offset=offset)),
                     1e-12*offset[idim]))
 
-        from meshmode.mesh.processing import glue_mesh_boundaries
-        periodic_mesh = glue_mesh_boundaries(mesh, glued_boundary_mappings)
+        periodic_mesh = glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols)
 
         return periodic_mesh
     else:
@@ -1350,9 +1351,11 @@ def transform(x):
         from meshmode.mesh.tools import AffineMap
         aff_map = AffineMap(matrix, center - matrix @ center)
 
-        from meshmode.mesh.processing import glue_mesh_boundaries
-        periodic_mesh = glue_mesh_boundaries(mesh,
-            glued_boundary_mappings=[("-theta", "+theta", aff_map, 1e-12)])
+        from meshmode.mesh.processing import (
+            glue_mesh_boundaries, BoundaryPairMapping)
+        periodic_mesh = glue_mesh_boundaries(
+            mesh, bdry_pair_mappings_and_tols=[
+                (BoundaryPairMapping("-theta", "+theta", aff_map), 1e-12)])
 
         return periodic_mesh
     else:
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 61d6c9e99..3ef22bc8c 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -24,6 +24,8 @@
 from numbers import Real
 from typing import Optional, Union
 
+from dataclasses import dataclass
+
 import numpy as np
 import numpy.linalg as la
 
@@ -35,8 +37,12 @@
     InterPartitionAdjacencyGroup
 )
 
+from meshmode.mesh.tools import AffineMap
+
 
 __doc__ = """
+.. autoclass:: BoundaryPairMapping
+
 .. autofunction:: find_group_indices
 .. autofunction:: partition_mesh
 .. autofunction:: find_volume_mesh_element_orientations
@@ -895,7 +901,6 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
 # something like this using numpy constructs
 def _match_vertices(
         mesh, src_vertex_indices, tgt_vertex_indices, aff_map=None, tol=1e-12):
-    from meshmode.mesh.tools import AffineMap
     if aff_map is None:
         aff_map = AffineMap()
 
@@ -941,6 +946,35 @@ def _match_vertices(
 
 # {{{ boundary face matching
 
+@dataclass(frozen=True)
+class BoundaryPairMapping:
+    """
+    Represents an affine mapping from one boundary to another.
+
+    .. attribute:: from_btag
+
+        The tag of one boundary.
+
+    .. attribute:: to_btag
+
+        The tag of the other boundary.
+
+    .. attribute:: aff_map
+
+        An :class:`meshmode.AffineMap` that maps points on boundary *from_btag* into
+        points on boundary *to_btag*.
+    """
+    from_btag: int
+    to_btag: int
+    aff_map: AffineMap
+
+    def inverted(self):
+        return BoundaryPairMapping(
+            self.to_btag,
+            self.from_btag,
+            self.aff_map.inverted())
+
+
 def _get_boundary_face_ids(mesh, btag):
     from meshmode.mesh import _FaceIDs
     face_ids_per_boundary_group = []
@@ -983,7 +1017,10 @@ def _get_face_vertex_indices(mesh, face_ids):
     return np.stack(face_vertex_indices_per_group)
 
 
-def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
+def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
+    btag_m = bdry_pair_mapping.from_btag
+    btag_n = bdry_pair_mapping.to_btag
+
     bdry_m_face_ids = _get_boundary_face_ids(mesh, btag_m)
     bdry_n_face_ids = _get_boundary_face_ids(mesh, btag_n)
 
@@ -1006,14 +1043,14 @@ def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
 
     matched_bdry_n_vertex_indices = _match_vertices(
         mesh, bdry_m_vertex_indices, bdry_n_vertex_indices,
-        aff_map=aff_map, tol=tol)
+        aff_map=bdry_pair_mapping.aff_map, tol=tol)
 
     unmatched_bdry_m_vertex_indices = bdry_m_vertex_indices[
         np.where(matched_bdry_n_vertex_indices < 0)[0]]
     nunmatched = len(unmatched_bdry_m_vertex_indices)
     if nunmatched > 0:
         vertices = mesh.vertices[:, unmatched_bdry_m_vertex_indices]
-        mapped_vertices = aff_map(vertices)
+        mapped_vertices = bdry_pair_mapping.aff_map(vertices)
         raise RuntimeError(
             f"unable to match vertices between boundaries {btag_m} and {btag_n}.\n"
             + "Unmatched vertices (original -> mapped):\n"
@@ -1054,22 +1091,21 @@ def _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol):
 
 # {{{ boundary gluing
 
-def _complete_glued_boundary_mappings(partial_glued_boundary_mappings):
+def _complete_boundary_pairs(partial_bdry_pair_mappings_and_tols):
     partial_btag_pairs = {
-        (btag_m, btag_n)
-        for btag_m, btag_n, _, _ in partial_glued_boundary_mappings}
+        (mapping.from_btag, mapping.to_btag)
+        for mapping, _ in partial_bdry_pair_mappings_and_tols}
 
-    glued_boundary_mappings = []
-    for btag_m, btag_n, aff_map, tol in partial_glued_boundary_mappings:
-        glued_boundary_mappings.append((btag_m, btag_n, aff_map, tol))
-        if (btag_n, btag_m) not in partial_btag_pairs:
-            inv_aff_map = aff_map.inverted()
-            glued_boundary_mappings.append((btag_n, btag_m, inv_aff_map, tol))
+    bdry_pair_mappings_and_tols = []
+    for mapping, tol in partial_bdry_pair_mappings_and_tols:
+        bdry_pair_mappings_and_tols.append((mapping, tol))
+        if (mapping.to_btag, mapping.from_btag) not in partial_btag_pairs:
+            bdry_pair_mappings_and_tols.append((mapping.inverted(), tol))
 
-    return glued_boundary_mappings
+    return bdry_pair_mappings_and_tols
 
 
-def glue_mesh_boundaries(mesh, glued_boundary_mappings):
+def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
     """
     Create a new mesh from *mesh* in which one or more pairs of boundaries are
     "glued" together such that the boundary surfaces become part of the interior
@@ -1080,25 +1116,24 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
     operates only on facial adjacency; any existing nodal adjacency in *mesh* is
     ignored/invalidated.
 
-    :arg glued_boundary_mappings: a :class:`list` of tuples
-        *(btag_m, btag_n, aff_map, tol)* which each specify a mapping between two
-        boundaries in *mesh* that should be glued together. *aff_map* is a
-        :class:`~meshmode.AffineMap` that represents the affine mapping from the
-        vertices of boundary *btag_m* into the vertices of boundary *btag_n*. *tol*
-        is the tolerance allowed between the vertex coordinates of *btag_n* and the
-        transformed vertex coordinates of *btag_m* when attempting to match the two.
+    :arg bdry_pair_mappings_and_tols: a :class:`list` of tuples *(mapping, tol)*,
+        where *mapping* is a :class:`BoundaryPairMapping` instance that specifies
+        a mapping between two boundaries in *mesh* that should be glued together,
+        and *tol* is the allowed tolerance between the transformed vertex
+        coordinates of the first boundary and the vertex coordinates of the second
+        boundary when attempting to match the two.
     """
-    glued_boundary_mappings = _complete_glued_boundary_mappings(
-        glued_boundary_mappings)
+    bdry_pair_mappings_and_tols = _complete_boundary_pairs(
+        bdry_pair_mappings_and_tols)
 
     glued_btags = {
         btag
-        for btag_m, btag_n, _, _ in glued_boundary_mappings
-        for btag in (btag_m, btag_n)}
+        for mapping, _ in bdry_pair_mappings_and_tols
+        for btag in (mapping.from_btag, mapping.to_btag)}
 
     face_id_pairs_for_mapping = [
-        _match_boundary_faces(mesh, btag_m, btag_n, aff_map, tol)
-        for btag_m, btag_n, aff_map, tol in glued_boundary_mappings]
+        _match_boundary_faces(mesh, mapping, tol)
+        for mapping, tol in bdry_pair_mappings_and_tols]
 
     from meshmode.mesh import InteriorAdjacencyGroup, BoundaryAdjacencyGroup
 
@@ -1110,7 +1145,7 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
             if not isinstance(fagrp, BoundaryAdjacencyGroup)
             or fagrp.boundary_tag not in glued_btags]
 
-        for imap, (_, _, aff_map, _) in enumerate(glued_boundary_mappings):
+        for imap, (mapping, _) in enumerate(bdry_pair_mappings_and_tols):
             mapping_face_id_pairs = face_id_pairs_for_mapping[imap]
             belongs_to_group = mapping_face_id_pairs[0].groups == igrp
             for ineighbor_grp in range(len(mesh.groups)):
@@ -1129,7 +1164,7 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings):
                         element_faces=element_faces,
                         neighbors=neighbors,
                         neighbor_faces=neighbor_faces,
-                        aff_map=aff_map))
+                        aff_map=mapping.aff_map))
 
         facial_adjacency_groups.append(fagrp_list)
 
@@ -1205,7 +1240,6 @@ def affine_map(mesh,
     if b is not None and b.shape != (mesh.ambient_dim,):
         raise ValueError(f"b has shape '{b.shape}' for a {mesh.ambient_dim}d mesh")
 
-    from meshmode.mesh.tools import AffineMap
     f = AffineMap(A, b)
 
     vertices = f(mesh.vertices)
diff --git a/test/test_mesh.py b/test/test_mesh.py
index 672f7801b..9db3ebba2 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -902,10 +902,11 @@ def test_glued_mesh():
     map_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]), center)
     map_upper_to_lower = _get_rotation(-np.pi/2, np.array([0, 0, 1]), center)
 
-    from meshmode.mesh.processing import glue_mesh_boundaries
-    mesh = glue_mesh_boundaries(orig_mesh,
-        glued_boundary_mappings=[
-            ("-theta", "+theta", map_lower_to_upper, 1e-12)
+    from meshmode.mesh.processing import (
+        glue_mesh_boundaries, BoundaryPairMapping)
+    mesh = glue_mesh_boundaries(
+        orig_mesh, bdry_pair_mappings_and_tols=[
+            (BoundaryPairMapping("-theta", "+theta", map_lower_to_upper), 1e-12)
         ])
 
     int_grps = [
@@ -970,10 +971,14 @@ def test_glued_mesh_matrix_only():
     matrix_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]))
     matrix_upper_to_lower = _get_rotation(-np.pi/2, np.array([0, 0, 1]))
 
-    from meshmode.mesh.processing import glue_mesh_boundaries
-    mesh = glue_mesh_boundaries(orig_mesh,
-        glued_boundary_mappings=[
-            ("-theta", "+theta", AffineMap(matrix=matrix_lower_to_upper), 1e-12)
+    map_lower_to_upper = AffineMap(matrix=matrix_lower_to_upper)
+    map_upper_to_lower = AffineMap(matrix=matrix_upper_to_lower)
+
+    from meshmode.mesh.processing import (
+        glue_mesh_boundaries, BoundaryPairMapping)
+    mesh = glue_mesh_boundaries(
+        orig_mesh, bdry_pair_mappings_and_tols=[
+            (BoundaryPairMapping("-theta", "+theta", map_lower_to_upper), 1e-12)
         ])
 
     int_grps = [
@@ -983,8 +988,8 @@ def test_glued_mesh_matrix_only():
     lower_grp = int_grps[1]
     upper_grp = int_grps[2]
 
-    assert lower_grp.aff_map == AffineMap(matrix=matrix_lower_to_upper)
-    assert upper_grp.aff_map == AffineMap(matrix=matrix_upper_to_lower)
+    assert lower_grp.aff_map == map_lower_to_upper
+    assert upper_grp.aff_map == map_upper_to_lower
 
 
 def test_glued_mesh_offset_only():
@@ -994,10 +999,14 @@ def test_glued_mesh_offset_only():
     offset_lower_to_upper = np.array([0, 0, 1])
     offset_upper_to_lower = np.array([0, 0, -1])
 
-    from meshmode.mesh.processing import glue_mesh_boundaries
-    mesh = glue_mesh_boundaries(orig_mesh,
-        glued_boundary_mappings=[
-            ("-z", "+z", AffineMap(offset=offset_lower_to_upper), 1e-12)
+    map_lower_to_upper = AffineMap(offset=offset_lower_to_upper)
+    map_upper_to_lower = AffineMap(offset=offset_upper_to_lower)
+
+    from meshmode.mesh.processing import (
+        glue_mesh_boundaries, BoundaryPairMapping)
+    mesh = glue_mesh_boundaries(
+        orig_mesh, bdry_pair_mappings_and_tols=[
+            (BoundaryPairMapping("-z", "+z", map_lower_to_upper), 1e-12)
         ])
 
     int_grps = [
@@ -1007,8 +1016,8 @@ def test_glued_mesh_offset_only():
     lower_grp = int_grps[1]
     upper_grp = int_grps[2]
 
-    assert lower_grp.aff_map == AffineMap(offset=offset_lower_to_upper)
-    assert upper_grp.aff_map == AffineMap(offset=offset_upper_to_lower)
+    assert lower_grp.aff_map == map_lower_to_upper
+    assert upper_grp.aff_map == map_upper_to_lower
 
 # }}}
 

From bb0a2eaf7ec44e4e88b8f245e75bf43d940f798e Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Wed, 6 Oct 2021 13:53:13 -0500
Subject: [PATCH 14/20] split mapping_face_id_pairs into two variables

---
 meshmode/mesh/processing.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 3ef22bc8c..210adbb09 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1146,17 +1146,17 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
             or fagrp.boundary_tag not in glued_btags]
 
         for imap, (mapping, _) in enumerate(bdry_pair_mappings_and_tols):
-            mapping_face_id_pairs = face_id_pairs_for_mapping[imap]
-            belongs_to_group = mapping_face_id_pairs[0].groups == igrp
+            face_ids, neighbor_face_ids = face_id_pairs_for_mapping[imap]
+            belongs_to_group = face_ids.groups == igrp
             for ineighbor_grp in range(len(mesh.groups)):
                 indices, = np.where(
                     belongs_to_group
-                    & (mapping_face_id_pairs[1].groups == ineighbor_grp))
+                    & (neighbor_face_ids.groups == ineighbor_grp))
                 if len(indices) > 0:
-                    elements = mapping_face_id_pairs[0].elements[indices]
-                    element_faces = mapping_face_id_pairs[0].faces[indices]
-                    neighbors = mapping_face_id_pairs[1].elements[indices]
-                    neighbor_faces = mapping_face_id_pairs[1].faces[indices]
+                    elements = face_ids.elements[indices]
+                    element_faces = face_ids.faces[indices]
+                    neighbors = neighbor_face_ids.elements[indices]
+                    neighbor_faces = neighbor_face_ids.faces[indices]
                     fagrp_list.append(InteriorAdjacencyGroup(
                         igroup=igrp,
                         ineighbor_group=ineighbor_grp,

From b640be2a5fa4319bec205db9bec31a3b3b2ca8e8 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Wed, 6 Oct 2021 15:03:35 -0500
Subject: [PATCH 15/20] remove sorting matching faces by group

not needed anymore, I think
---
 meshmode/mesh/processing.py | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 210adbb09..ad2bacfb3 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1073,18 +1073,16 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
 
     assert face_index_pairs.shape[1] == nfaces
 
-    order = np.argsort(face_index_pairs[0, :])
-
     from meshmode.mesh import _FaceIDs
     return (
         _FaceIDs(
-            groups=face_ids.groups[face_index_pairs[0, order]],
-            elements=face_ids.elements[face_index_pairs[0, order]],
-            faces=face_ids.faces[face_index_pairs[0, order]]),
+            groups=face_ids.groups[face_index_pairs[0, :]],
+            elements=face_ids.elements[face_index_pairs[0, :]],
+            faces=face_ids.faces[face_index_pairs[0, :]]),
         _FaceIDs(
-            groups=face_ids.groups[face_index_pairs[1, order]],
-            elements=face_ids.elements[face_index_pairs[1, order]],
-            faces=face_ids.faces[face_index_pairs[1, order]]))
+            groups=face_ids.groups[face_index_pairs[1, :]],
+            elements=face_ids.elements[face_index_pairs[1, :]],
+            faces=face_ids.faces[face_index_pairs[1, :]]))
 
 # }}}
 

From 70bfaf028708ddd4fdafcecd08d32f11b4558d5b Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Thu, 7 Oct 2021 14:41:06 -0500
Subject: [PATCH 16/20] specify order of intra-match indices in
 _match_faces_by_vertices

---
 meshmode/mesh/__init__.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 759b282a9..6b07fdbad 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1159,7 +1159,9 @@ def _match_faces_by_vertices(groups, face_ids, vertex_index_map_func=None):
         indices. Must accept multidimensional arrays as input and return an array
         of the same shape.
     :returns: A :class:`numpy.ndarray` of shape ``(2, nmatches)`` of indices into
-        *face_ids*.
+        *face_ids*. The ordering of the matches returned is unspecified. For a given
+        match, however, the first index will correspond to the face that occurs first
+        in *face_ids*.
     """
     if vertex_index_map_func is None:
         def vertex_index_map_func(vertices):

From d4e9ab4d68d3a08077a76723d22c624b633acf92 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Thu, 7 Oct 2021 15:19:29 -0500
Subject: [PATCH 17/20] document _match_boundary_faces

---
 meshmode/mesh/processing.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index ad2bacfb3..fa1da3d19 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1018,6 +1018,23 @@ def _get_face_vertex_indices(mesh, face_ids):
 
 
 def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
+    """
+    Given a :class:`BoundaryPairMapping` *bdry_pair_mapping*, return the
+    correspondence between faces of the two boundaries (expressed as a pair of
+    :class:`meshmode.mesh._FaceIDs`).
+
+    :arg mesh: The mesh containing the boundaries.
+    :arg bdry_pair_mapping: A :class:`BoundaryPairMapping` specifying the boundaries
+        whose faces are to be matched.
+    :arg tol: The allowed tolerance between the transformed vertex coordinates of
+        the first boundary and the vertex coordinates of the second boundary.
+    :returns: A pair of :class:`meshmode.mesh._FaceIDs`, each having a number of
+        entries equal to the number of faces in the boundary, that represents the
+        correspondence between the two boundaries' faces. The first element in the
+        pair contains faces from boundary *bdry_pair_mapping.from_btag*, and the
+        second contains faces from boundary *bdry_pair_mapping.to_btag*. The order
+        of the faces is unspecified.
+    """
     btag_m = bdry_pair_mapping.from_btag
     btag_n = bdry_pair_mapping.to_btag
 
@@ -1073,6 +1090,9 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
 
     assert face_index_pairs.shape[1] == nfaces
 
+    # Since the first boundary's faces come before the second boundary's in
+    # face_ids, the first boundary's faces should all be in the first row of the
+    # result of _match_faces_by_vertices
     from meshmode.mesh import _FaceIDs
     return (
         _FaceIDs(

From 2ddc9ec9b441788d27a1df0bbf217091f7004235 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Thu, 7 Oct 2021 16:16:39 -0500
Subject: [PATCH 18/20] don't do face matching twice

---
 meshmode/mesh/processing.py | 77 +++++++++++++++++++++----------------
 1 file changed, 43 insertions(+), 34 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index fa1da3d19..ee76d3c14 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -968,12 +968,6 @@ class BoundaryPairMapping:
     to_btag: int
     aff_map: AffineMap
 
-    def inverted(self):
-        return BoundaryPairMapping(
-            self.to_btag,
-            self.from_btag,
-            self.aff_map.inverted())
-
 
 def _get_boundary_face_ids(mesh, btag):
     from meshmode.mesh import _FaceIDs
@@ -1109,20 +1103,6 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
 
 # {{{ boundary gluing
 
-def _complete_boundary_pairs(partial_bdry_pair_mappings_and_tols):
-    partial_btag_pairs = {
-        (mapping.from_btag, mapping.to_btag)
-        for mapping, _ in partial_bdry_pair_mappings_and_tols}
-
-    bdry_pair_mappings_and_tols = []
-    for mapping, tol in partial_bdry_pair_mappings_and_tols:
-        bdry_pair_mappings_and_tols.append((mapping, tol))
-        if (mapping.to_btag, mapping.from_btag) not in partial_btag_pairs:
-            bdry_pair_mappings_and_tols.append((mapping.inverted(), tol))
-
-    return bdry_pair_mappings_and_tols
-
-
 def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
     """
     Create a new mesh from *mesh* in which one or more pairs of boundaries are
@@ -1139,16 +1119,28 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
         a mapping between two boundaries in *mesh* that should be glued together,
         and *tol* is the allowed tolerance between the transformed vertex
         coordinates of the first boundary and the vertex coordinates of the second
-        boundary when attempting to match the two.
+        boundary when attempting to match the two. Pass at most one mapping for each
+        unique (order-independent) pair of boundaries.
     """
-    bdry_pair_mappings_and_tols = _complete_boundary_pairs(
-        bdry_pair_mappings_and_tols)
-
     glued_btags = {
         btag
         for mapping, _ in bdry_pair_mappings_and_tols
         for btag in (mapping.from_btag, mapping.to_btag)}
 
+    btag_to_index = {btag: i for i, btag in enumerate(glued_btags)}
+
+    glued_btag_pairs = set()
+    for mapping, _ in bdry_pair_mappings_and_tols:
+        if btag_to_index[mapping.from_btag] < btag_to_index[mapping.to_btag]:
+            btag_pair = (mapping.from_btag, mapping.to_btag)
+        else:
+            btag_pair = (mapping.to_btag, mapping.from_btag)
+        if btag_pair in glued_btag_pairs:
+            raise ValueError(
+                "multiple mappings detected for boundaries "
+                f"{btag_pair[0]} and {btag_pair[1]}.")
+        glued_btag_pairs.add(btag_pair)
+
     face_id_pairs_for_mapping = [
         _match_boundary_faces(mesh, mapping, tol)
         for mapping, tol in bdry_pair_mappings_and_tols]
@@ -1164,17 +1156,21 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
             or fagrp.boundary_tag not in glued_btags]
 
         for imap, (mapping, _) in enumerate(bdry_pair_mappings_and_tols):
-            face_ids, neighbor_face_ids = face_id_pairs_for_mapping[imap]
-            belongs_to_group = face_ids.groups == igrp
+            bdry_m_face_ids, bdry_n_face_ids = face_id_pairs_for_mapping[imap]
+            bdry_m_belongs_to_group = bdry_m_face_ids.groups == igrp
+            bdry_n_belongs_to_group = bdry_n_face_ids.groups == igrp
             for ineighbor_grp in range(len(mesh.groups)):
-                indices, = np.where(
-                    belongs_to_group
-                    & (neighbor_face_ids.groups == ineighbor_grp))
-                if len(indices) > 0:
-                    elements = face_ids.elements[indices]
-                    element_faces = face_ids.faces[indices]
-                    neighbors = neighbor_face_ids.elements[indices]
-                    neighbor_faces = neighbor_face_ids.faces[indices]
+                bdry_m_indices, = np.where(
+                    bdry_m_belongs_to_group
+                    & (bdry_n_face_ids.groups == ineighbor_grp))
+                bdry_n_indices, = np.where(
+                    bdry_n_belongs_to_group
+                    & (bdry_m_face_ids.groups == ineighbor_grp))
+                if len(bdry_m_indices) > 0:
+                    elements = bdry_m_face_ids.elements[bdry_m_indices]
+                    element_faces = bdry_m_face_ids.faces[bdry_m_indices]
+                    neighbors = bdry_n_face_ids.elements[bdry_m_indices]
+                    neighbor_faces = bdry_n_face_ids.faces[bdry_m_indices]
                     fagrp_list.append(InteriorAdjacencyGroup(
                         igroup=igrp,
                         ineighbor_group=ineighbor_grp,
@@ -1183,6 +1179,19 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
                         neighbors=neighbors,
                         neighbor_faces=neighbor_faces,
                         aff_map=mapping.aff_map))
+                if len(bdry_n_indices) > 0:
+                    elements = bdry_n_face_ids.elements[bdry_n_indices]
+                    element_faces = bdry_n_face_ids.faces[bdry_n_indices]
+                    neighbors = bdry_m_face_ids.elements[bdry_n_indices]
+                    neighbor_faces = bdry_m_face_ids.faces[bdry_n_indices]
+                    fagrp_list.append(InteriorAdjacencyGroup(
+                        igroup=igrp,
+                        ineighbor_group=ineighbor_grp,
+                        elements=elements,
+                        element_faces=element_faces,
+                        neighbors=neighbors,
+                        neighbor_faces=neighbor_faces,
+                        aff_map=mapping.aff_map.inverted()))
 
         facial_adjacency_groups.append(fagrp_list)
 

From 06565f4d3f242236580430e5efde4e17031b85d6 Mon Sep 17 00:00:00 2001
From: Matthew Smith <mjsmith6@illinois.edu>
Date: Fri, 8 Oct 2021 14:17:22 -0500
Subject: [PATCH 19/20] add numpy vertex matching for small cases

---
 meshmode/mesh/processing.py | 93 +++++++++++++++++++++++--------------
 test/test_mesh.py           |  5 +-
 2 files changed, 61 insertions(+), 37 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index ee76d3c14..d5a95db64 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -897,47 +897,66 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
 
 # {{{ vertex matching
 
-# FIXME: This tree-based approach is probably slow; see if there's a way to do
-# something like this using numpy constructs
 def _match_vertices(
-        mesh, src_vertex_indices, tgt_vertex_indices, aff_map=None, tol=1e-12):
+        mesh, src_vertex_indices, tgt_vertex_indices, *, aff_map=None, tol=1e-12,
+        use_tree=None):
     if aff_map is None:
         aff_map = AffineMap()
 
+    if use_tree is None:
+        # Rough empirical guess for when the tree version becomes faster
+        use_tree = len(tgt_vertex_indices) >= 2**13
+
     src_vertices = mesh.vertices[:, src_vertex_indices]
     tgt_vertices = mesh.vertices[:, tgt_vertex_indices]
 
-    tgt_vertex_bboxes = np.stack((
-        tgt_vertices - tol,
-        tgt_vertices + tol))
+    mapped_src_vertices = aff_map(src_vertices)
 
-    from pytools.spatial_btree import SpatialBinaryTreeBucket
-    tree = SpatialBinaryTreeBucket(
-        np.min(tgt_vertex_bboxes[0], axis=1),
-        np.max(tgt_vertex_bboxes[1], axis=1))
-    for ivertex in range(len(tgt_vertex_indices)):
-        tree.insert(ivertex, tgt_vertex_bboxes[:, :, ivertex])
+    if use_tree:
+        tgt_vertex_bboxes = np.stack((
+            tgt_vertices - tol,
+            tgt_vertices + tol))
+
+        from pytools.spatial_btree import SpatialBinaryTreeBucket
+        tree = SpatialBinaryTreeBucket(
+            np.min(tgt_vertex_bboxes[0], axis=1),
+            np.max(tgt_vertex_bboxes[1], axis=1))
+        for ivertex in range(len(tgt_vertex_indices)):
+            tree.insert(ivertex, tgt_vertex_bboxes[:, :, ivertex])
+
+        matched_tgt_vertices = np.full(len(src_vertex_indices), -1)
+        for ivertex in range(len(src_vertex_indices)):
+            mapped_src_vertex = mapped_src_vertices[:, ivertex]
+            matches = np.array(list(tree.generate_matches(mapped_src_vertex)))
+            match_bboxes = tgt_vertex_bboxes[:, :, matches]
+            in_bbox = np.all(
+                (mapped_src_vertex[:, np.newaxis] >= match_bboxes[0, :, :])
+                & (mapped_src_vertex[:, np.newaxis] <= match_bboxes[1, :, :]),
+                axis=0)
+            candidate_indices = matches[in_bbox]
+            if len(candidate_indices) == 0:
+                continue
+            displacements = (
+                mapped_src_vertex.reshape(-1, 1)
+                - tgt_vertices[:, candidate_indices])
+            distances_sq = np.sum(displacements**2, axis=0)
+            matched_tgt_vertices[ivertex] = (
+                tgt_vertex_indices[candidate_indices[np.argmin(distances_sq)]])
 
-    mapped_src_vertices = aff_map(src_vertices)
+    else:
+        displacements = (
+            mapped_src_vertices.reshape(mesh.dim, -1, 1)
+            - tgt_vertices.reshape(mesh.dim, 1, -1))
+        distances_sq = np.sum(displacements**2, axis=0)
 
-    matched_tgt_vertices = np.full(len(src_vertex_indices), -1)
-    for ivertex in range(len(src_vertex_indices)):
-        mapped_src_vertex = mapped_src_vertices[:, ivertex]
-        matches = np.array(list(tree.generate_matches(mapped_src_vertex)))
-        match_bboxes = tgt_vertex_bboxes[:, :, matches]
-        in_bbox = np.all(
-            (mapped_src_vertex[:, np.newaxis] >= match_bboxes[0, :, :])
-            & (mapped_src_vertex[:, np.newaxis] <= match_bboxes[1, :, :]),
-            axis=0)
-        candidate_indices = matches[in_bbox]
-        if len(candidate_indices) == 0:
-            continue
-        displacement = (
-            mapped_src_vertex.reshape(-1, 1)
-            - tgt_vertices[:, candidate_indices])
-        distance_sq = np.sum(displacement**2, axis=0)
-        matched_tgt_vertices[ivertex] = (
-            tgt_vertex_indices[candidate_indices[np.argmin(distance_sq)]])
+        vertex_indices, = np.indices((len(src_vertex_indices),))
+        min_distance_sq_indices = np.argmin(distances_sq, axis=1)
+        min_distances_sq = distances_sq[vertex_indices, min_distance_sq_indices]
+
+        matched_tgt_vertices = np.where(
+            min_distances_sq < tol**2,
+            tgt_vertex_indices[min_distance_sq_indices],
+            -1)
 
     return matched_tgt_vertices
 
@@ -1011,7 +1030,7 @@ def _get_face_vertex_indices(mesh, face_ids):
     return np.stack(face_vertex_indices_per_group)
 
 
-def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
+def _match_boundary_faces(mesh, bdry_pair_mapping, tol, *, use_tree=None):
     """
     Given a :class:`BoundaryPairMapping` *bdry_pair_mapping*, return the
     correspondence between faces of the two boundaries (expressed as a pair of
@@ -1022,6 +1041,8 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
         whose faces are to be matched.
     :arg tol: The allowed tolerance between the transformed vertex coordinates of
         the first boundary and the vertex coordinates of the second boundary.
+    :arg use_tree: Optional argument indicating whether to use a spatial binary
+        search tree or a (quadratic) numpy algorithm when matching vertices.
     :returns: A pair of :class:`meshmode.mesh._FaceIDs`, each having a number of
         entries equal to the number of faces in the boundary, that represents the
         correspondence between the two boundaries' faces. The first element in the
@@ -1054,7 +1075,7 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
 
     matched_bdry_n_vertex_indices = _match_vertices(
         mesh, bdry_m_vertex_indices, bdry_n_vertex_indices,
-        aff_map=bdry_pair_mapping.aff_map, tol=tol)
+        aff_map=bdry_pair_mapping.aff_map, tol=tol, use_tree=use_tree)
 
     unmatched_bdry_m_vertex_indices = bdry_m_vertex_indices[
         np.where(matched_bdry_n_vertex_indices < 0)[0]]
@@ -1103,7 +1124,7 @@ def _match_boundary_faces(mesh, bdry_pair_mapping, tol):
 
 # {{{ boundary gluing
 
-def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
+def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols, *, use_tree=None):
     """
     Create a new mesh from *mesh* in which one or more pairs of boundaries are
     "glued" together such that the boundary surfaces become part of the interior
@@ -1121,6 +1142,8 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
         coordinates of the first boundary and the vertex coordinates of the second
         boundary when attempting to match the two. Pass at most one mapping for each
         unique (order-independent) pair of boundaries.
+    :arg use_tree: Optional argument indicating whether to use a spatial binary
+        search tree or a (quadratic) numpy algorithm when matching vertices.
     """
     glued_btags = {
         btag
@@ -1142,7 +1165,7 @@ def glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols):
         glued_btag_pairs.add(btag_pair)
 
     face_id_pairs_for_mapping = [
-        _match_boundary_faces(mesh, mapping, tol)
+        _match_boundary_faces(mesh, mapping, tol, use_tree=use_tree)
         for mapping, tol in bdry_pair_mappings_and_tols]
 
     from meshmode.mesh import InteriorAdjacencyGroup, BoundaryAdjacencyGroup
diff --git a/test/test_mesh.py b/test/test_mesh.py
index 9db3ebba2..e1320d96c 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -893,7 +893,8 @@ def test_cube_icosphere(actx_factory, order, visualize=True):
 
 # {{{ mesh boundary gluing
 
-def test_glued_mesh():
+@pytest.mark.parametrize("use_tree", [False, True])
+def test_glued_mesh(use_tree):
     n = 4
     center = (1, 2, 3)
 
@@ -907,7 +908,7 @@ def test_glued_mesh():
     mesh = glue_mesh_boundaries(
         orig_mesh, bdry_pair_mappings_and_tols=[
             (BoundaryPairMapping("-theta", "+theta", map_lower_to_upper), 1e-12)
-        ])
+        ], use_tree=use_tree)
 
     int_grps = [
         fagrp for fagrp in mesh.facial_adjacency_groups[0]

From e28f594afbecd86f808df2ae7f880c2230167a68 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= <inform@tiker.net>
Date: Wed, 20 Oct 2021 12:27:03 -0500
Subject: [PATCH 20/20] Bump down cutoff for tree algorithm in face matching

To avoid giant temporaries.
---
 meshmode/mesh/processing.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index d5a95db64..577aad504 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -904,8 +904,11 @@ def _match_vertices(
         aff_map = AffineMap()
 
     if use_tree is None:
-        # Rough empirical guess for when the tree version becomes faster
-        use_tree = len(tgt_vertex_indices) >= 2**13
+        # Empirically, the tree version becomes faster at 2**13.
+        # The temporary (displacements) below at that size requires
+        # 1.6GB, which seems like a lot. Capping at 2**11 instead,
+        # which requires a more reasonable 100M.
+        use_tree = len(tgt_vertex_indices) >= 2**11
 
     src_vertices = mesh.vertices[:, src_vertex_indices]
     tgt_vertices = mesh.vertices[:, tgt_vertex_indices]