Add mesh boundary gluing (aka periodic)

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)
 

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,

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())
 

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)