diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py
index b3f2d3a07..12c263d69 100644
--- a/meshmode/discretization/connection/opposite_face.py
+++ b/meshmode/discretization/connection/opposite_face.py
@@ -42,7 +42,8 @@ def thaw_to_numpy(actx, array):
 def _make_cross_face_batches(actx,
         tgt_bdry_discr, src_bdry_discr,
         i_tgt_grp, i_src_grp,
-        tgt_bdry_element_indices, src_bdry_element_indices):
+        tgt_bdry_element_indices, src_bdry_element_indices,
+        tgt_aff_transforms=None, src_aff_transforms=None):
 
     if tgt_bdry_discr.dim == 0:
         return [InterpolationBatch(
@@ -52,15 +53,25 @@ def _make_cross_face_batches(actx,
             result_unit_nodes=src_bdry_discr.groups[i_src_grp].unit_nodes,
             to_element_face=None)]
 
-    tgt_bdry_nodes = np.array([
+    def transform(aff_transforms, x):
+        mats, vecs = aff_transforms if aff_transforms is not None else (None, None)
+        if mats is not None:
+            result = np.einsum("edi,ien->den", mats, x)
+        else:
+            result = x
+        if vecs is not None:
+            result = result + vecs.T.reshape(vecs.shape[1], -1, 1)
+        return result
+
+    tgt_bdry_nodes = transform(tgt_aff_transforms, np.array([
         thaw_to_numpy(actx, ary[i_tgt_grp])[tgt_bdry_element_indices]
         for ary in tgt_bdry_discr.nodes(cached=False)
-        ])
+        ]))
 
-    src_bdry_nodes = np.array([
+    src_bdry_nodes = transform(src_aff_transforms, np.array([
         thaw_to_numpy(actx, ary[i_src_grp])[src_bdry_element_indices]
         for ary in src_bdry_discr.nodes(cached=False)
-        ])
+        ]))
 
     tol = 1e4 * np.finfo(tgt_bdry_nodes.dtype).eps
 
@@ -449,6 +460,19 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn):
                         vbc_tgt_grp_face_batch.to_element_indices
                         )[vbc_used_els]
 
+                mats = (
+                    adj.aff_transform_mats[adj_tgt_flags, :, :]
+                    if adj.aff_transform_mats is not None
+                    else None)
+                vecs = (
+                    adj.aff_transform_vecs[adj_tgt_flags, :]
+                    if adj.aff_transform_vecs is not None
+                    else None)
+                tgt_aff_transforms = (
+                    (mats, vecs)
+                    if mats is not None or vecs is not None
+                    else None)
+
                 # find src_bdry_element_indices
 
                 src_vol_element_indices = adj.neighbors[adj_tgt_flags]
@@ -488,7 +512,8 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn):
                         bdry_discr, bdry_discr,
                         i_tgt_grp, i_src_grp,
                         tgt_bdry_element_indices,
-                        src_bdry_element_indices)
+                        src_bdry_element_indices,
+                        tgt_aff_transforms=tgt_aff_transforms)
                 groups[i_tgt_grp].extend(batches)
 
     from meshmode.discretization.connection import (
@@ -561,6 +586,15 @@ def make_partition_connection(actx, *, local_bdry_conn, i_local_part,
         i_local_vol_elems = rem_ipag.partition_neighbors[indices]
         i_local_faces = rem_ipag.neighbor_faces[indices]
 
+        remote_aff_mats = (
+            rem_ipag.aff_transform_mats[indices, :, :]
+            if rem_ipag.aff_transform_mats is not None
+            else None)
+        remote_aff_vecs = (
+            rem_ipag.aff_transform_vecs[indices, :]
+            if rem_ipag.aff_transform_vecs is not None
+            else None)
+
         del indices
 
         i_local_grps = find_group_indices(local_vol_groups, i_local_vol_elems)
@@ -603,13 +637,27 @@ def make_partition_connection(actx, *, local_bdry_conn, i_local_part,
                     i_remote_faces[local_indices]]
             assert (matched_remote_bdry_el_indices >= 0).all()
 
+            mats = (
+                remote_aff_mats[local_indices, :, :]
+                if remote_aff_mats is not None
+                else None)
+            vecs = (
+                remote_aff_vecs[local_indices, :]
+                if remote_aff_vecs is not None
+                else None)
+            src_aff_transforms = (
+                (mats, vecs)
+                if mats is not None or vecs is not None
+                else None)
+
             # }}}
 
             grp_batches = _make_cross_face_batches(actx,
                         local_bdry_conn.to_discr, remote_bdry_discr,
                         i_local_grp, rem_ipag.igroup,
                         matched_local_bdry_el_indices,
-                        matched_remote_bdry_el_indices)
+                        matched_remote_bdry_el_indices,
+                        src_aff_transforms=src_aff_transforms)
 
             part_batches[i_local_grp].extend(grp_batches)
 
diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py
index bbc9e762f..60dc994aa 100644
--- a/meshmode/interop/firedrake/mesh.py
+++ b/meshmode/interop/firedrake/mesh.py
@@ -386,7 +386,7 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology,
                                              new_ext_neighbor_faces))
         ext_neighbors = np.concatenate((ext_neighbors, new_ext_neighbors))
 
-    exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor=None,
+    exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor_group=None,
                                         elements=ext_elements,
                                         element_faces=ext_element_faces,
                                         neighbors=ext_neighbors,
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 583536863..4bbba72b6 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -468,10 +468,35 @@ class FacialAdjacencyGroup(Record):
 
         Zero if ``neighbors[i]`` is negative.
 
+    .. attribute:: aff_transform_mats
+
+        ``np.float64 [nfagrp_elements, ambient_dim, ambient_dim]`` or None.
+        ``aff_transform_mats[iface, :, :]`` gives the matrix part of the affine
+        mapping from face *iface* to the corresponding neighbor face.
+
+    .. attribute:: aff_transform_vecs
+
+        ``np.float64 [nfagrp_elements, ambient_dim]`` or None.
+        ``aff_transform_vecs[iface, :]`` gives the vector part of the affine
+        mapping from face *iface* to the corresponding neighbor face.
+
     .. automethod:: __eq__
     .. automethod:: __ne__
     """
 
+    def __init__(self, igroup, ineighbor_group, *,
+            elements, element_faces,
+            neighbors, neighbor_faces,
+            aff_transform_mats=None, aff_transform_vecs=None,
+            **kwargs):
+        Record.__init__(self,
+            igroup=igroup, ineighbor_group=ineighbor_group,
+            elements=elements, element_faces=element_faces,
+            neighbors=neighbors, neighbor_faces=neighbor_faces,
+            aff_transform_mats=aff_transform_mats,
+            aff_transform_vecs=aff_transform_vecs,
+            **kwargs)
+
     def __eq__(self, other):
         return (
                 type(self) == type(other)
@@ -481,6 +506,8 @@ def __eq__(self, other):
                 and np.array_equal(self.element_faces, other.element_faces)
                 and np.array_equal(self.neighbors, other.neighbors)
                 and np.array_equal(self.neighbor_faces, other.neighbor_faces)
+                and np.array_equal(self.aff_transform_mats, other.aff_transform_mats)
+                and np.array_equal(self.aff_transform_vecs, other.aff_transform_vecs)
                 )
 
     def __ne__(self, other):
@@ -554,9 +581,37 @@ class InterPartitionAdjacencyGroup(FacialAdjacencyGroup):
         If ``neighbor_partitions[i]`` is negative, ``elements[i]`` is on a true
         boundary and is not connected to any other :class:``Mesh``.
 
+    .. attribute:: aff_transform_mats
+
+        ``np.float64 [nfagrp_elements, ambient_dim, ambient_dim]`` or None.
+        ``aff_transform_mats[iface, :, :]`` gives the matrix part of the affine
+        mapping from face *iface* to the corresponding neighbor face.
+
+    .. attribute:: aff_transform_vecs
+
+        ``np.float64 [nfagrp_elements, ambient_dim]`` or None.
+        ``aff_transform_vecs[iface, :]`` gives the vector part of the affine
+        mapping from face *iface* to the corresponding neighbor face.
+
     .. versionadded:: 2017.1
     """
 
+    def __init__(self, igroup, ineighbor_group, *,
+            elements, element_faces,
+            neighbors, neighbor_faces,
+            partition_neighbors, neighbor_partitions,
+            aff_transform_mats=None, aff_transform_vecs=None,
+            **kwargs):
+        FacialAdjacencyGroup.__init__(self,
+            igroup=igroup, ineighbor_group=ineighbor_group,
+            elements=elements, element_faces=element_faces,
+            neighbors=neighbors, neighbor_faces=neighbor_faces,
+            aff_transform_mats=aff_transform_mats,
+            aff_transform_vecs=aff_transform_vecs,
+            partition_neighbors=partition_neighbors,
+            neighbor_partitions=neighbor_partitions,
+            **kwargs)
+
     def __eq__(self, other):
         return (super.__eq__(self, other)
             and np.array_equal(self.partition_neighbors, other.partition_neighbors)
@@ -1302,6 +1357,14 @@ def fagrp_params_str(fagrp):
                     "element_faces": _numpy_array_as_python(fagrp.element_faces),
                     "neighbors": _numpy_array_as_python(fagrp.neighbors),
                     "neighbor_faces": _numpy_array_as_python(fagrp.neighbor_faces),
+                    "aff_transform_mats": (
+                        _numpy_array_as_python(fagrp.aff_transform_mats)
+                        if fagrp.aff_transform_mats is not None
+                        else None),
+                    "aff_transform_vecs": (
+                        _numpy_array_as_python(fagrp.aff_transform_vecs)
+                        if fagrp.aff_transform_vecs is not None
+                        else None),
                     }
             return ",\n    ".join(f"{k}={v}" for k, v in params.items())
 
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index aa3d78f75..49ec03099 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -69,6 +69,8 @@
 .. autofunction:: generate_box_mesh
 .. autofunction:: generate_regular_rect_mesh
 .. autofunction:: generate_warped_rect_mesh
+.. autofunction:: generate_annular_cylinder_slice_mesh
+.. autofunction:: generate_periodic_annular_cylinder_slice_mesh
 
 Tools for Iterative Refinement
 ------------------------------
@@ -1213,6 +1215,68 @@ def m(x):
 # }}}
 
 
+# {{{ generate_annular_cylinder_slice_mesh
+
+def generate_annular_cylinder_slice_mesh(n, center, inner_radius, outer_radius):
+    r"""
+    Generate a slice of a 3D annular cylinder for
+    :math:`\theta \in [-\frac{\pi}{4}, \frac{\pi}{4}]`.
+    """
+    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)))
+
+    return mesh
+
+# }}}
+
+
+# {{{ generate_periodic_annular_cylinder_slice_mesh
+
+def generate_periodic_annular_cylinder_slice_mesh(
+        n, center, inner_radius, outer_radius):
+    r"""
+    Generate a slice of a 3D annular cylinder for
+    :math:`\theta \in [-\frac{\pi}{4}, \frac{\pi}{4}]` that is periodic in
+    :math:`\theta`.
+    """
+    nonperiodic_mesh = generate_annular_cylinder_slice_mesh(
+        n, center, inner_radius, outer_radius)
+
+    from meshmode.mesh.processing import _get_rotation_matrix_from_angle_and_axis
+    mat = _get_rotation_matrix_from_angle_and_axis(np.pi/2, np.array([0, 0, 1]))
+    vec = center - mat @ center
+
+    from meshmode.mesh.processing import glue_mesh_boundaries
+    mesh = glue_mesh_boundaries(nonperiodic_mesh,
+        glued_boundary_mappings=[("-theta", "+theta", (mat, vec), 1e-12)])
+
+    return mesh
+
+# }}}
+
+
 # {{{ warp_and_refine_until_resolved
 
 @log_process(logger)
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 695250a75..c7f7da4d0 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -40,6 +40,7 @@
 .. autofunction:: find_bounding_box
 .. autofunction:: merge_disjoint_meshes
 .. autofunction:: split_mesh_groups
+.. autofunction:: glue_mesh_boundaries
 
 .. autofunction:: map_mesh
 .. autofunction:: affine_map
@@ -65,6 +66,14 @@ def find_group_indices(groups, meshwide_elems):
     return grps
 
 
+def _repeat_matrix(mat, n):
+    """
+    Return an array of shape ``(n, mat.shape[0], mat.shape[1])`` consisting of *n*
+    copies of *mat*.
+    """
+    return np.repeat(mat.reshape(1, mat.shape[0], mat.shape[1]), n, axis=0)
+
+
 # {{{ partition_mesh
 
 def _compute_global_elem_to_part_elem(part_per_element, parts, element_id_dtype):
@@ -224,6 +233,14 @@ def _create_local_to_local_adjacency_groups(mesh, global_elem_to_part_elem,
                 neighbors = global_elem_to_part_elem[facial_adj.neighbors[
                             adj_indices] + elem_base_j] - part_elem_base_j
                 neighbor_faces = facial_adj.neighbor_faces[adj_indices]
+                mats = (
+                    facial_adj.aff_transform_mats[adj_indices, :, :]
+                    if facial_adj.aff_transform_mats is not None
+                    else None)
+                vecs = (
+                    facial_adj.aff_transform_vecs[adj_indices, :]
+                    if facial_adj.aff_transform_vecs is not None
+                    else None)
 
                 from meshmode.mesh import FacialAdjacencyGroup
                 local_to_local_adjacency_groups[i_part_grp][j_part_grp] =\
@@ -232,7 +249,9 @@ def _create_local_to_local_adjacency_groups(mesh, global_elem_to_part_elem,
                                         elements=elements,
                                         element_faces=element_faces,
                                         neighbors=neighbors,
-                                        neighbor_faces=neighbor_faces)
+                                        neighbor_faces=neighbor_faces,
+                                        aff_transform_mats=mats,
+                                        aff_transform_vecs=vecs)
 
     return local_to_local_adjacency_groups
 
@@ -263,12 +282,22 @@ class _NonLocalAdjacencyData:
     .. attribute:: neighbor_faces
 
        The index of the shared face inside the remote element.
+
+    .. attribute:: aff_transform_mats
+
+       The matrix part of the affine mapping from local to remote element.
+
+    .. attribute:: aff_transform_vecs
+
+       The vector part of the affine mapping from local to remote element.
     """
     elements: np.ndarray
     element_faces: np.ndarray
     neighbor_parts: np.ndarray
     global_neighbors: np.ndarray
     neighbor_faces: np.ndarray
+    aff_transform_mats: np.ndarray
+    aff_transform_vecs: np.ndarray
 
 
 def _collect_nonlocal_adjacency_data(mesh, part_per_elem, global_elem_to_part_elem,
@@ -326,17 +355,51 @@ def _collect_nonlocal_adjacency_data(mesh, part_per_elem, global_elem_to_part_el
                 global_neighbors = facial_adj.neighbors[adj_indices] + elem_base_j
                 neighbor_parts = part_per_elem[global_neighbors]
                 neighbor_faces = facial_adj.neighbor_faces[adj_indices]
-
-                pairwise_adj.append(_NonLocalAdjacencyData(elements, element_faces,
-                            neighbor_parts, global_neighbors, neighbor_faces))
+                mats = (
+                    facial_adj.aff_transform_mats[adj_indices, :, :]
+                    if facial_adj.aff_transform_mats is not None
+                    else None)
+                vecs = (
+                    facial_adj.aff_transform_vecs[adj_indices, :]
+                    if facial_adj.aff_transform_vecs is not None
+                    else None)
+
+                pairwise_adj.append(
+                    _NonLocalAdjacencyData(
+                        elements, element_faces, neighbor_parts, global_neighbors,
+                        neighbor_faces, mats, vecs))
 
         if pairwise_adj:
+            if any([
+                    adj.aff_transform_mats is not None
+                    for adj in pairwise_adj]):
+                eye = np.eye(mesh.ambient_dim, dtype=np.float64)
+                concatenated_mats = np.concatenate([
+                    adj.aff_transform_mats
+                    if adj.aff_transform_mats is not None
+                    else _repeat_matrix(eye, len(adj.elements))
+                    for adj in pairwise_adj])
+            else:
+                concatenated_mats = None
+            if any([
+                    adj.aff_transform_vecs is not None
+                    for adj in pairwise_adj]):
+                concatenated_vecs = np.concatenate([
+                    adj.aff_transform_vecs
+                    if adj.aff_transform_vecs is not None
+                    else np.zeros((len(adj.elements), mesh.ambient_dim),
+                        dtype=np.float64)
+                    for adj in pairwise_adj])
+            else:
+                concatenated_vecs = None
             nonlocal_adj_data[i_part_grp] = _NonLocalAdjacencyData(
                 np.concatenate([adj.elements for adj in pairwise_adj]),
                 np.concatenate([adj.element_faces for adj in pairwise_adj]),
                 np.concatenate([adj.neighbor_parts for adj in pairwise_adj]),
                 np.concatenate([adj.global_neighbors for adj in pairwise_adj]),
-                np.concatenate([adj.neighbor_faces for adj in pairwise_adj]))
+                np.concatenate([adj.neighbor_faces for adj in pairwise_adj]),
+                concatenated_mats,
+                concatenated_vecs)
 
     return nonlocal_adj_data
 
@@ -452,6 +515,8 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
             neighbors = np.empty(0, dtype=mesh.element_id_dtype)
             neighbor_elements = np.empty(0, dtype=mesh.element_id_dtype)
             neighbor_faces = np.empty(0, dtype=mesh.face_id_dtype)
+            mats = None
+            vecs = None
 
         elif bdry is None:
             # Non-local adjacency only
@@ -466,6 +531,8 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
             neighbors = -flags
             neighbor_elements = global_elem_to_neighbor_elem[nl.global_neighbors]
             neighbor_faces = nl.neighbor_faces
+            mats = nl.aff_transform_mats
+            vecs = nl.aff_transform_vecs
 
         elif nl is None:
             # Boundary only
@@ -476,6 +543,8 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
             neighbors = bdry.neighbors
             neighbor_elements = np.full(nfaces, -1, dtype=mesh.element_id_dtype)
             neighbor_faces = np.zeros(nfaces, dtype=mesh.face_id_dtype)
+            mats = None
+            vecs = None
 
         else:
             # Both; need to merge together
@@ -488,6 +557,14 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
             neighbors = np.empty(nfaces, dtype=mesh.element_id_dtype)
             neighbor_elements = np.empty(nfaces, dtype=mesh.element_id_dtype)
             neighbor_faces = np.empty(nfaces, dtype=mesh.face_id_dtype)
+            mats = (
+                _repeat_matrix(np.eye(mesh.ambient_dim, dtype=np.float64), nfaces)
+                if nl.aff_transform_mats is not None
+                else None)
+            vecs = (
+                np.zeros((nfaces, mesh.ambient_dim), dtype=np.float64)
+                if nl.aff_transform_vecs is not None
+                else None)
 
             # Combine lists of elements/faces and sort to assist in merging
             combined_elements = np.concatenate((nl.elements, bdry.elements))
@@ -509,6 +586,10 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
             neighbor_elements[nonlocal_indices] = global_elem_to_neighbor_elem[
                         nl.global_neighbors]
             neighbor_faces[nonlocal_indices] = nl.neighbor_faces
+            if nl.aff_transform_mats is not None:
+                mats[nonlocal_indices] = nl.aff_transform_mats
+            if nl.aff_transform_vecs is not None:
+                vecs[nonlocal_indices] = nl.aff_transform_vecs
 
             # Merge boundary part
             bdry_indices = np.where(perm >= nnonlocal)[0]
@@ -525,7 +606,9 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element,
                     element_faces=element_faces, neighbors=neighbors,
                     neighbor_partitions=neighbor_parts,
                     partition_neighbors=neighbor_elements,
-                    neighbor_faces=neighbor_faces))
+                    neighbor_faces=neighbor_faces,
+                    aff_transform_mats=mats,
+                    aff_transform_vecs=vecs))
 
     return inter_partition_adj_groups
 
@@ -1019,12 +1102,460 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
 # }}}
 
 
+# {{{ mesh boundary gluing
+
+def _get_bdry_face_ids(mesh, btag):
+    btag_bit = mesh.boundary_tag_bit(btag)
+
+    from meshmode.mesh import _FaceIDs
+    face_ids_per_group = []
+    for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups):
+        bdry_grp = fagrp_map.get(None)
+        if bdry_grp is not None:
+            bdry_face_indices, = np.where((-bdry_grp.neighbors & btag_bit) != 0)
+            grp_face_ids = _FaceIDs(
+                groups=np.zeros(len(bdry_face_indices), dtype=int) + igrp,
+                elements=bdry_grp.elements[bdry_face_indices],
+                faces=bdry_grp.element_faces[bdry_face_indices])
+        else:
+            grp_face_ids = _FaceIDs(
+                groups=np.empty(0, dtype=int),
+                elements=np.empty(0, dtype=mesh.element_id_dtype),
+                faces=np.empty(0, dtype=mesh.face_id_dtype))
+        face_ids_per_group.append(grp_face_ids)
+
+    from meshmode.mesh import _concatenate_face_ids
+    return _concatenate_face_ids(face_ids_per_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):
+        is_grp = face_ids.groups == igrp
+        face_vertex_indices = np.full(
+            (np.count_nonzero(is_grp), max_face_vertices), -1,
+            dtype=mesh.vertex_id_dtype)
+        for fid, ref_fvi in enumerate(grp.face_vertex_indices()):
+            is_face = is_grp & (face_ids.faces == fid)
+            is_face_grp = face_ids.faces[is_grp] == fid
+            face_vertex_indices[is_face_grp, :len(ref_fvi)] = (
+                grp.vertex_indices[face_ids.elements[is_face], :][:, ref_fvi])
+        face_vertex_indices_per_group.append(face_vertex_indices)
+
+    return np.stack(face_vertex_indices_per_group)
+
+
+def _compute_face_indices_from_mask(mask):
+    # Order by face, then by element
+    indices = np.cumsum(mask).reshape(mask.shape) - 1
+    indices[~mask] = -1
+    return indices
+
+
+def _match_boundary_faces(mesh, glued_boundary_mappings):
+    face_id_pairs_for_mapping = []
+
+    for btag_m, btag_n, aff_transform, tol in glued_boundary_mappings:
+        bdry_m_face_ids = _get_bdry_face_ids(mesh, btag_m)
+        bdry_n_face_ids = _get_bdry_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]
+
+        nvertices = single_valued((
+            len(bdry_m_vertex_indices),
+            len(bdry_n_vertex_indices)))
+
+        bdry_m_vertices = mesh.vertices[:, bdry_m_vertex_indices]
+        bdry_n_vertices = mesh.vertices[:, bdry_n_vertex_indices]
+
+        # FIXME: This approach is probably slow; see if there's a way to do
+        # something like this using numpy constructs
+        bdry_n_bbox = (
+            np.min(bdry_n_vertices, axis=1),
+            np.max(bdry_n_vertices, axis=1))
+        from pytools.spatial_btree import SpatialBinaryTreeBucket
+        tree = SpatialBinaryTreeBucket(bdry_n_bbox[0], bdry_n_bbox[1])
+        bdry_n_vertex_bboxes = np.stack((
+            bdry_n_vertices - tol,
+            bdry_n_vertices + tol))
+        for ivertex in range(nvertices):
+            tree.insert(ivertex, bdry_n_vertex_bboxes[:, :, ivertex])
+
+        mat, vec = aff_transform
+        if mat is not None:
+            mapped_bdry_m_vertices = mat @ bdry_m_vertices
+        else:
+            mapped_bdry_m_vertices = bdry_m_vertices
+        if vec is not None:
+            mapped_bdry_m_vertices = mapped_bdry_m_vertices + vec.reshape(-1, 1)
+
+        equivalent_vertices = np.empty((2, nvertices), dtype=mesh.element_id_dtype)
+        for ivertex in range(nvertices):
+            bdry_m_index = bdry_m_vertex_indices[ivertex]
+            mapped_bdry_m_vertex = mapped_bdry_m_vertices[:, ivertex]
+            matches = np.array(list(tree.generate_matches(mapped_bdry_m_vertex)))
+            match_bboxes = bdry_n_vertex_bboxes[:, :, matches]
+            in_bbox = np.all(
+                (mapped_bdry_m_vertex[:, np.newaxis] >= match_bboxes[0, :, :])
+                & (mapped_bdry_m_vertex[:, np.newaxis] <= match_bboxes[1, :, :]),
+                axis=0)
+            candidate_indices = matches[in_bbox]
+            if len(candidate_indices) == 0:
+                raise RuntimeError("failed to find a matching vertex for vertex "
+                    f"{bdry_m_index} at {mapped_bdry_m_vertex}")
+            displacement = (
+                mapped_bdry_m_vertex.reshape(-1, 1)
+                - bdry_n_vertices[:, candidate_indices])
+            distance_sq = np.sum(displacement**2, axis=0)
+            bdry_n_index = bdry_n_vertex_indices[
+                candidate_indices[np.argmin(distance_sq)]]
+            equivalent_vertices[:, ivertex] = [bdry_m_index, bdry_n_index]
+
+        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[equivalent_vertices[0, :]] = equivalent_vertices[1, :]
+
+        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
+
+        from meshmode.mesh import _FaceIDs
+        order = np.argsort(face_index_pairs[0, :])
+        face_id_pairs_for_mapping.append((
+            _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]])))
+
+    return face_id_pairs_for_mapping
+
+
+def _translate_boundary_pair_adjacency_into_group_pair_adjacency(mesh,
+        glued_boundary_mappings, face_id_pairs_for_mapping):
+    face_id_pairs_for_group_pair = []
+    mapping_indices_for_group_pair = []
+
+    for igrp, grp in enumerate(mesh.groups):
+        face_id_pairs_grp_map = {}
+        mapping_indices_grp_map = {}
+
+        connected_groups = np.unique(np.concatenate([
+            face_id_pairs[1].groups[face_id_pairs[0].groups == igrp]
+            for face_id_pairs in face_id_pairs_for_mapping]))
+
+        for ineighbor_grp in connected_groups:
+            adj_data_for_mapping = []
+            face_has_neighbor = np.full((grp.nfaces, grp.nelements), False)
+
+            for imap in range(len(glued_boundary_mappings)):
+                mapping_face_id_pairs = face_id_pairs_for_mapping[imap]
+                is_grp_pair = (
+                    (mapping_face_id_pairs[0].groups == igrp)
+                    & (mapping_face_id_pairs[1].groups == ineighbor_grp))
+                elements = mapping_face_id_pairs[0].elements[is_grp_pair]
+                element_faces = mapping_face_id_pairs[0].faces[is_grp_pair]
+                neighbors = mapping_face_id_pairs[1].elements[is_grp_pair]
+                neighbor_faces = mapping_face_id_pairs[1].faces[is_grp_pair]
+                adj_data_for_mapping.append(
+                    (elements, element_faces, neighbors, neighbor_faces))
+                face_has_neighbor[element_faces, elements] = True
+
+            translated_indices = _compute_face_indices_from_mask(face_has_neighbor)
+            nfaces = np.max(translated_indices) + 1
+
+            from meshmode.mesh import _FaceIDs
+            face_id_pairs = (
+                _FaceIDs(
+                    groups=np.zeros(nfaces, dtype=int) + igrp,
+                    elements=np.empty(nfaces, dtype=mesh.element_id_dtype),
+                    faces=np.empty(nfaces, dtype=mesh.face_id_dtype)),
+                _FaceIDs(
+                    groups=np.zeros(nfaces, dtype=int) + ineighbor_grp,
+                    elements=np.empty(nfaces, dtype=mesh.element_id_dtype),
+                    faces=np.empty(nfaces, dtype=mesh.face_id_dtype)))
+
+            mapping_indices = np.empty(nfaces, dtype=int)
+
+            for imap, (elements, element_faces, neighbors, neighbor_faces) in (
+                    enumerate(adj_data_for_mapping)):
+                indices = translated_indices[element_faces, elements]
+                face_id_pairs[0].elements[indices] = elements
+                face_id_pairs[0].faces[indices] = element_faces
+                face_id_pairs[1].elements[indices] = neighbors
+                face_id_pairs[1].faces[indices] = neighbor_faces
+                mapping_indices[indices] = imap
+
+            face_id_pairs_grp_map[ineighbor_grp] = face_id_pairs
+            mapping_indices_grp_map[ineighbor_grp] = mapping_indices
+
+        face_id_pairs_for_group_pair.append(face_id_pairs_grp_map)
+        mapping_indices_for_group_pair.append(mapping_indices_grp_map)
+
+    return face_id_pairs_for_group_pair, mapping_indices_for_group_pair
+
+
+def _construct_glued_mesh(mesh, glued_boundary_mappings,
+        face_id_pairs_for_group_pair, mapping_indices_for_group_pair):
+    glued_btags = (
+        set(btag_m for btag_m, _, _, _ in glued_boundary_mappings)
+        | set(btag_n for _, btag_n, _, _ in glued_boundary_mappings))
+
+    boundary_tags = [
+        btag for btag in mesh.boundary_tags
+        if btag not in glued_btags]
+
+    btag_to_index = {btag: i for i, btag in enumerate(boundary_tags)}
+
+    def boundary_tag_bit(btag):
+        from meshmode.mesh import _boundary_tag_bit
+        return _boundary_tag_bit(boundary_tags, btag_to_index, btag)
+
+    mapping_has_mat = np.array([
+        mat is not None for _, _, (mat, _), _ in glued_boundary_mappings])
+    mapping_has_vec = np.array([
+        vec is not None for _, _, (_, vec), _ in glued_boundary_mappings])
+
+    mats_for_mapping = np.stack([
+        mat if mat is not None
+        else np.eye(mesh.ambient_dim, dtype=np.float64)
+        for _, _, (mat, _), _ in glued_boundary_mappings])
+    vecs_for_mapping = np.stack([
+        vec if vec is not None
+        else np.zeros(mesh.ambient_dim, dtype=np.float64)
+        for _, _, (_, vec), _ in glued_boundary_mappings])
+
+    from meshmode.mesh import FacialAdjacencyGroup
+
+    facial_adjacency_groups = []
+
+    for igrp, grp in enumerate(mesh.groups):
+        fagrp_map = {}
+
+        old_fagrp_map = mesh.facial_adjacency_groups[igrp]
+        face_id_pairs_grp_map = face_id_pairs_for_group_pair[igrp]
+        mapping_indices_grp_map = mapping_indices_for_group_pair[igrp]
+
+        connected_groups = (
+            set(
+                ineighbor_grp for ineighbor_grp in old_fagrp_map.keys()
+                if ineighbor_grp is not None)
+            | set(ineighbor_grp for ineighbor_grp in face_id_pairs_grp_map.keys()))
+
+        for ineighbor_grp in connected_groups:
+            face_has_neighbor = np.full((grp.nfaces, grp.nelements), False)
+            old_adj = old_fagrp_map.get(ineighbor_grp)
+            if old_adj is not None:
+                face_has_neighbor[old_adj.element_faces, old_adj.elements] = True
+            grp_pair_face_ids = face_id_pairs_grp_map.get(ineighbor_grp)
+            if grp_pair_face_ids is not None:
+                face_ids = grp_pair_face_ids[0]
+                face_has_neighbor[face_ids.faces, face_ids.elements] = True
+
+            merged_indices = _compute_face_indices_from_mask(face_has_neighbor)
+            nfaces = np.max(merged_indices) + 1
+
+            mapping_indices = mapping_indices_grp_map[ineighbor_grp]
+
+            elements = np.empty(nfaces, dtype=mesh.element_id_dtype)
+            element_faces = np.empty(nfaces, dtype=mesh.face_id_dtype)
+            neighbors = np.empty(nfaces, dtype=mesh.element_id_dtype)
+            neighbor_faces = np.empty(nfaces, dtype=mesh.face_id_dtype)
+
+            mats = (
+                _repeat_matrix(np.eye(mesh.ambient_dim, dtype=np.float64), nfaces)
+                if (
+                    old_adj.aff_transform_mats is not None
+                    or np.any(mapping_has_mat[mapping_indices]))
+                else None)
+            vecs = (
+                np.zeros((nfaces, mesh.ambient_dim), dtype=np.float64)
+                if (
+                    old_adj.aff_transform_vecs is not None
+                    or np.any(mapping_has_vec[mapping_indices]))
+                else None)
+
+            if old_adj is not None:
+                indices = merged_indices[old_adj.element_faces, old_adj.elements]
+                elements[indices] = old_adj.elements
+                element_faces[indices] = old_adj.element_faces
+                neighbors[indices] = old_adj.neighbors
+                neighbor_faces[indices] = old_adj.neighbor_faces
+                if old_adj.aff_transform_mats is not None:
+                    mats[indices, :, :] = old_adj.aff_transform_mats
+                if old_adj.aff_transform_vecs is not None:
+                    vecs[indices, :] = old_adj.aff_transform_vecs
+
+            if grp_pair_face_ids is not None:
+                face_ids = grp_pair_face_ids[0]
+                neighbor_face_ids = grp_pair_face_ids[1]
+                indices = merged_indices[face_ids.faces, face_ids.elements]
+                elements[indices] = face_ids.elements
+                element_faces[indices] = face_ids.faces
+                neighbors[indices] = neighbor_face_ids.elements
+                neighbor_faces[indices] = neighbor_face_ids.faces
+                if mats is not None:
+                    mats[indices, :, :] = mats_for_mapping[mapping_indices, :, :]
+                if vecs is not None:
+                    vecs[indices, :] = vecs_for_mapping[mapping_indices, :]
+
+            fagrp_map[ineighbor_grp] = FacialAdjacencyGroup(
+                igroup=igrp,
+                ineighbor_group=ineighbor_grp,
+                elements=elements,
+                element_faces=element_faces,
+                neighbors=neighbors,
+                neighbor_faces=neighbor_faces,
+                aff_transform_mats=mats,
+                aff_transform_vecs=vecs)
+
+        is_bdry = np.full((grp.nfaces, grp.nelements), True)
+        for ineighbor_grp in connected_groups:
+            adj = fagrp_map.get(ineighbor_grp)
+            is_bdry[adj.element_faces, adj.elements] = False
+
+        if np.any(is_bdry):
+            old_bdry_grp = old_fagrp_map[None]
+            indices, = np.where(is_bdry[
+                old_bdry_grp.element_faces,
+                old_bdry_grp.elements])
+
+            old_flags = -old_bdry_grp.neighbors[indices]
+            flags = old_flags.copy()
+            for btag in mesh.boundary_tags:
+                old_btag_bit = mesh.boundary_tag_bit(btag)
+                flags &= ~old_btag_bit
+            for btag in boundary_tags:
+                old_btag_bit = mesh.boundary_tag_bit(btag)
+                btag_bit = boundary_tag_bit(btag)
+                flags[(old_flags & old_btag_bit) != 0] |= btag_bit
+
+            fagrp_map[None] = FacialAdjacencyGroup(
+                igroup=igrp,
+                ineighbor_group=None,
+                elements=old_bdry_grp.elements[indices],
+                element_faces=old_bdry_grp.element_faces[indices],
+                neighbors=-flags,
+                neighbor_faces=old_bdry_grp.neighbor_faces[indices])
+
+        facial_adjacency_groups.append(fagrp_map)
+
+    return mesh.copy(
+        boundary_tags=boundary_tags,
+        facial_adjacency_groups=facial_adjacency_groups)
+
+
+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_transform, tol)* which each specify a mapping between
+        two boundaries in *mesh* that should be glued together. *aff_transform* is a
+        tuple *(mat, vec)* 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.
+    """
+    mapped_btags = set(
+        (btag_m, btag_n)
+        for btag_m, btag_n, _, _ in glued_boundary_mappings)
+
+    glued_boundary_mappings_both_ways = []
+
+    for btag_m, btag_n, aff_transform, tol in glued_boundary_mappings:
+        transform_mat = (
+            np.array(aff_transform[0])
+            if aff_transform is not None and aff_transform[0] is not None
+            else None)
+        transform_vec = (
+            np.array(aff_transform[1])
+            if aff_transform is not None and aff_transform[1] is not None
+            else None)
+        glued_boundary_mappings_both_ways.append(
+            (btag_m, btag_n, (transform_mat, transform_vec), tol))
+        if (btag_n, btag_m) not in mapped_btags:
+            if transform_mat is not None:
+                inv_transform_mat = la.inv(transform_mat)
+                inv_transform_vec = (
+                    -inv_transform_mat @ transform_vec
+                    if transform_vec is not None
+                    else None)
+            else:
+                inv_transform_mat = None
+                inv_transform_vec = (
+                    # https://github.com/PyCQA/pylint/issues/4608
+                    -transform_vec  # pylint: disable=invalid-unary-operand-type
+                    if transform_vec is not None
+                    else None)
+            glued_boundary_mappings_both_ways.append(
+                (btag_n, btag_m, (inv_transform_mat, inv_transform_vec), tol))
+
+    face_id_pairs_for_mapping = _match_boundary_faces(mesh,
+        glued_boundary_mappings_both_ways)
+
+    face_id_pairs_for_group_pair, mapping_indices_for_group_pair = (
+        _translate_boundary_pair_adjacency_into_group_pair_adjacency(
+            mesh, glued_boundary_mappings_both_ways, face_id_pairs_for_mapping))
+
+    glued_mesh = _construct_glued_mesh(mesh, glued_boundary_mappings_both_ways,
+        face_id_pairs_for_group_pair, mapping_indices_for_group_pair)
+
+    return glued_mesh
+
+# }}}
+
+
 # {{{ map
 
 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_transforms = any([
+            adj.aff_transform_mats is not None or adj.aff_transform_vecs is not None
+            for fagrp_map in mesh.facial_adjacency_groups
+            for adj in fagrp_map.values()])
+        if has_adj_transforms:
+            raise ValueError("cannot apply a general map to a mesh that has "
+                "transforms 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")
@@ -1072,20 +1603,79 @@ 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")
 
-    if b is not None:
-        b = b.reshape(-1, 1)
-
     def f(x):
         z = x
         if A is not None:
             z = A @ z
 
         if b is not None:
-            z = z + b
+            z = z + b.reshape(-1, 1)
 
         return z
 
-    return map_mesh(mesh, f)
+    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 adjancency 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
+        has_adj_transform_vecs = any([
+            adj.aff_transform_vecs is not None
+            for fagrp_map in mesh.facial_adjacency_groups
+            for adj in fagrp_map.values()])
+        if has_adj_transform_vecs:
+            facial_adjacency_groups = []
+            for old_fagrp_map in mesh.facial_adjacency_groups:
+                fagrp_map = {}
+                for ineighbor_grp, old_adj in old_fagrp_map.items():
+                    if ineighbor_grp is None:
+                        continue
+                    if old_adj.aff_transform_vecs is not None:
+                        if A is not None:
+                            aff_transform_vecs = (
+                                np.einsum("ij,fj->fi", A,
+                                    old_adj.aff_transform_vecs))
+                        else:
+                            aff_transform_vecs = old_adj.aff_transform_vecs.copy()
+                        if b is not None:
+                            mats = old_adj.aff_transform_mats
+                            aff_transform_vecs += b.reshape(1, -1) - mats @ b
+                        fagrp_map[ineighbor_grp] = old_adj.copy(
+                            aff_transform_vecs=aff_transform_vecs)
+                    else:
+                        fagrp_map[ineighbor_grp] = old_adj.copy()
+                facial_adjacency_groups.append(fagrp_map)
+        else:
+            facial_adjacency_groups = mesh.facial_adjacency_groups.copy()
+    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 56f331815..421da302d 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -200,6 +200,109 @@ 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_transforms(visualize=False):
+    orig_mesh = mgen.generate_periodic_annular_cylinder_slice_mesh(
+        4, (1, 2, 0), 0.5, 1)
+
+    if visualize:
+        from meshmode.mesh.visualization import write_vertex_vtk_file
+        write_vertex_vtk_file(orig_mesh, "affine_map_facial_adj_original.vtu")
+
+    def get_rotation(amount, axis, center=None):
+        from meshmode.mesh.processing import _get_rotation_matrix_from_angle_and_axis
+        mat = _get_rotation_matrix_from_angle_and_axis(amount, axis)
+        if center is None:
+            return mat
+        else:
+            # x0 + R @ (x - x0) = R @ x + (I - R) @ x0
+            vec = (np.eye(orig_mesh.ambient_dim) - mat) @ center
+            return mat, vec
+
+    mat_lower_to_upper = get_rotation(np.pi/2, np.array([0, 0, 1]))
+    mat_upper_to_lower = get_rotation(-np.pi/2, np.array([0, 0, 1]))
+
+    def find_matching_matrices(mats, query_mat, tol):
+        diff = np.abs(mats - query_mat[np.newaxis, :, :])
+        return np.where(np.einsum("fij->f", diff) < tol)[0]
+
+    tol = 1e-12
+
+    orig_adj = orig_mesh.facial_adjacency_groups[0][0]
+
+    periodic_lower_face_indices = find_matching_matrices(
+        orig_adj.aff_transform_mats, mat_lower_to_upper, tol)
+    periodic_upper_face_indices = find_matching_matrices(
+        orig_adj.aff_transform_mats, mat_upper_to_lower, tol)
+    nonperiodic_face_indices = find_matching_matrices(
+        orig_adj.aff_transform_mats, np.eye(orig_mesh.ambient_dim), tol)
+
+    assert len(periodic_lower_face_indices) > 0
+    assert len(periodic_upper_face_indices) > 0
+    assert len(nonperiodic_face_indices) > 0
+
+    from meshmode.mesh.processing import affine_map
+
+    def check_adj_aff_transforms(adj, *,
+            lower_to_upper_transform,
+            upper_to_lower_transform,
+            max_err=1e-12):
+        assert la.norm(
+            adj.aff_transform_mats[periodic_lower_face_indices, :, :]
+            - lower_to_upper_transform[0]) < max_err
+        assert la.norm(
+            adj.aff_transform_vecs[periodic_lower_face_indices, :]
+            - lower_to_upper_transform[1]) < max_err
+
+        assert la.norm(
+            adj.aff_transform_mats[periodic_upper_face_indices, :, :]
+            - upper_to_lower_transform[0]) < max_err
+        assert la.norm(
+            adj.aff_transform_vecs[periodic_upper_face_indices, :]
+            - upper_to_lower_transform[1]) < max_err
+
+        assert la.norm(
+            adj.aff_transform_mats[nonperiodic_face_indices, :, :]
+            - np.eye(orig_mesh.ambient_dim)) < max_err
+        assert la.norm(adj.aff_transform_vecs[nonperiodic_face_indices, :]) < max_err
+
+    # Matrix only
+    mesh = affine_map(orig_mesh, A=get_rotation(np.pi/2, np.array([0, 0, 1])))
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_matrix.vtu")
+
+    check_adj_aff_transforms(mesh.facial_adjacency_groups[0][0],
+        lower_to_upper_transform=get_rotation(
+            np.pi/2, np.array([0, 0, 1]), np.array([-2, 1, 0])),
+        upper_to_lower_transform=get_rotation(
+            -np.pi/2, np.array([0, 0, 1]), np.array([-2, 1, 0])))
+
+    # Vector only
+    mesh = affine_map(orig_mesh, b=np.array([0, -2, 0]))
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_vector.vtu")
+
+    check_adj_aff_transforms(mesh.facial_adjacency_groups[0][0],
+        lower_to_upper_transform=get_rotation(
+            np.pi/2, np.array([0, 0, 1]), np.array([1, 0, 0])),
+        upper_to_lower_transform=get_rotation(
+            -np.pi/2, np.array([0, 0, 1]), np.array([1, 0, 0])))
+
+    # Matrix and vector
+    transform = get_rotation(np.pi/2, np.array([0, 0, 1]), np.array([1, 1, 0]))
+    mesh = affine_map(orig_mesh, A=transform[0], b=transform[1])
+
+    if visualize:
+        write_vertex_vtk_file(mesh, "affine_map_facial_adj_matrix_and_vector.vtu")
+
+    check_adj_aff_transforms(mesh.facial_adjacency_groups[0][0],
+        lower_to_upper_transform=get_rotation(
+            np.pi/2, np.array([0, 0, 1]), np.array([0, 1, 0])),
+        upper_to_lower_transform=get_rotation(
+            -np.pi/2, np.array([0, 0, 1]), np.array([0, 1, 0])))
+
+
 @pytest.mark.parametrize("ambient_dim", [2, 3])
 def test_mesh_rotation(ambient_dim, visualize=False):
     order = 3
@@ -762,6 +865,158 @@ def test_quad_mesh_3d(mesh_name, order=3, visualize=False):
 # }}}
 
 
+# {{{ mesh boundary gluing
+
+def _get_rotation(amount, axis, center=None):
+    from meshmode.mesh.processing import _get_rotation_matrix_from_angle_and_axis
+    mat = _get_rotation_matrix_from_angle_and_axis(amount, axis)
+    if center is None:
+        return mat
+    else:
+        # x0 + R @ (x - x0) = R @ x + (I - R) @ x0
+        vec = (np.eye(3) - mat) @ center
+        return mat, vec
+
+
+def test_glued_mesh():
+    n = 4
+    center = (1, 2, 3)
+
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(n, center, 0.5, 1)
+
+    transform_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]), center)
+    transform_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", transform_lower_to_upper, 1e-12)
+        ])
+
+    assert "-theta" not in mesh.boundary_tags
+    assert "+theta" not in mesh.boundary_tags
+
+    def find_matching_transforms(mats, vecs, query_transform, tol):
+        mat_diff = np.abs(mats - query_transform[0][np.newaxis, :, :])
+        vec_diff = np.abs(vecs - query_transform[1][np.newaxis, :])
+        return np.where(
+            (np.einsum("fij->f", mat_diff) < tol)
+            & (np.einsum("fi->f", vec_diff) < tol))[0]
+
+    tol = 1e-12
+
+    orig_adj = orig_mesh.facial_adjacency_groups[0][0]
+    adj = mesh.facial_adjacency_groups[0][0]
+    bdry = mesh.facial_adjacency_groups[0][None]
+
+    assert adj.aff_transform_mats is not None
+    assert adj.aff_transform_vecs is not None
+
+    periodic_lower_face_indices = find_matching_transforms(
+        adj.aff_transform_mats, adj.aff_transform_vecs, transform_lower_to_upper,
+        tol)
+    periodic_upper_face_indices = find_matching_transforms(
+        adj.aff_transform_mats, adj.aff_transform_vecs, transform_upper_to_lower,
+        tol)
+    nonperiodic_face_indices = find_matching_transforms(
+        adj.aff_transform_mats, adj.aff_transform_vecs,
+        (np.eye(orig_mesh.ambient_dim), np.array([0, 0, 0])), tol)
+
+    assert len(periodic_lower_face_indices) == 2*n**2
+    assert len(periodic_upper_face_indices) == 2*n**2
+    assert len(nonperiodic_face_indices) == len(orig_adj.elements)
+    assert len(bdry.elements) == 8*n**2
+
+    occurrence_counts = np.zeros(
+        (mesh.groups[0].nfaces, mesh.nelements), dtype=int)
+
+    # Check that every face present in the facial adjacency has a match
+    occurrence_counts[adj.element_faces, adj.elements] += 1
+    occurrence_counts[adj.neighbor_faces, adj.neighbors] += 1
+    assert np.count_nonzero(occurrence_counts == 2) == len(adj.elements)
+
+    # Check that boundary faces and interior faces are disjoint
+    occurrence_counts[bdry.element_faces, bdry.elements] += 1
+    assert np.count_nonzero(occurrence_counts == 1) == len(bdry.elements)
+
+
+def test_partial_glued_mesh():
+    n = 4
+    orig_mesh = mgen.generate_annular_cylinder_slice_mesh(n, (0, 0, 0), 0.5, 1)
+
+    # Matrix only
+    mat_lower_to_upper = _get_rotation(np.pi/2, np.array([0, 0, 1]))
+    mat_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", (mat_lower_to_upper, None), 1e-12)
+        ])
+
+    orig_adj = orig_mesh.facial_adjacency_groups[0][0]
+    adj = mesh.facial_adjacency_groups[0][0]
+    bdry = mesh.facial_adjacency_groups[0][None]
+
+    assert adj.aff_transform_mats is not None
+    assert adj.aff_transform_vecs is None
+
+    def find_matching_matrices(mats, query_mat, tol):
+        diff = np.abs(mats - query_mat[np.newaxis, :, :])
+        return np.where(np.einsum("fij->f", diff) < tol)[0]
+
+    tol = 1e-12
+
+    periodic_lower_face_indices = find_matching_matrices(
+        adj.aff_transform_mats, mat_lower_to_upper, tol)
+    periodic_upper_face_indices = find_matching_matrices(
+        adj.aff_transform_mats, mat_upper_to_lower, tol)
+    nonperiodic_face_indices = find_matching_matrices(
+        adj.aff_transform_mats, np.eye(orig_mesh.ambient_dim), tol)
+
+    assert len(periodic_lower_face_indices) == 2*n**2
+    assert len(periodic_upper_face_indices) == 2*n**2
+    assert len(nonperiodic_face_indices) == len(orig_adj.elements)
+    assert len(bdry.elements) == 8*n**2
+
+    # Vector only
+    vec_lower_to_upper = np.array([0, 0, 1])
+    vec_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", (None, vec_lower_to_upper), 1e-12)
+        ])
+
+    orig_adj = orig_mesh.facial_adjacency_groups[0][0]
+    adj = mesh.facial_adjacency_groups[0][0]
+    bdry = mesh.facial_adjacency_groups[0][None]
+
+    assert adj.aff_transform_mats is None
+    assert adj.aff_transform_vecs is not None
+
+    def find_matching_vectors(vecs, query_vec, tol):
+        diff = np.abs(vecs - query_vec[np.newaxis, :])
+        return np.where(np.einsum("fi->f", diff) < tol)[0]
+
+    tol = 1e-12
+
+    periodic_lower_face_indices = find_matching_vectors(
+        adj.aff_transform_vecs, vec_lower_to_upper, tol)
+    periodic_upper_face_indices = find_matching_vectors(
+        adj.aff_transform_vecs, vec_upper_to_lower, tol)
+    nonperiodic_face_indices = find_matching_vectors(
+        adj.aff_transform_vecs, np.array([0, 0, 0]), tol)
+
+    assert len(periodic_lower_face_indices) == 2*n**2
+    assert len(periodic_upper_face_indices) == 2*n**2
+    assert len(nonperiodic_face_indices) == len(orig_adj.elements)
+    assert len(bdry.elements) == 8*n**2
+
+# }}}
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1:
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 6471eb7aa..a8d56c0db 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_periodic_annular_cylinder_slice_mesh(
+                mesh_par, (1, 2, 3), 0.5, 1)
+
             h = 1/mesh_par
         else:
             raise ValueError("mesh_name not recognized")