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): diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index dabb95c50..faf7a27b9 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 ------------------------------ @@ -880,14 +881,17 @@ def generate_urchin( # {{{ generate_box_mesh @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, +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. :arg 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. + :arg periodic: an optional tuple of :class:`bool` indicating whether + 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. @@ -946,6 +950,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 @@ -1078,12 +1085,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 = { @@ -1140,10 +1155,33 @@ 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.processing import ( + glue_mesh_boundaries, BoundaryPairMapping) + + from meshmode import AffineMap + bdry_pair_mappings_and_tols = [] + for idim in range(dim): + if periodic[idim]: + offset = np.zeros(dim, dtype=np.float64) + offset[idim] = axis_coords[idim][-1] - axis_coords[idim][0] + bdry_pair_mappings_and_tols.append(( + BoundaryPairMapping( + "periodic_-" + axes[idim], + "periodic_+" + axes[idim], + AffineMap(offset=offset)), + 1e-12*offset[idim])) + + periodic_mesh = glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols) + + return periodic_mesh + else: + return mesh + # }}} @@ -1152,6 +1190,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, @@ -1166,6 +1205,9 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), *, nelements_per_axis=None, number of elements along each axis. :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. 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`. @@ -1205,6 +1247,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) @@ -1267,6 +1310,60 @@ 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, BoundaryPairMapping) + periodic_mesh = glue_mesh_boundaries( + mesh, bdry_pair_mappings_and_tols=[ + (BoundaryPairMapping("-theta", "+theta", aff_map), 1e-12)]) + + return periodic_mesh + else: + return mesh + +# }}} + + # {{{ warp_and_refine_until_resolved @log_process(logger) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 9f00475c8..577aad504 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 @@ -45,6 +51,7 @@ .. autofunction:: find_bounding_box .. autofunction:: merge_disjoint_meshes .. autofunction:: split_mesh_groups +.. autofunction:: glue_mesh_boundaries .. autofunction:: map_mesh .. autofunction:: affine_map @@ -888,12 +895,357 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): # }}} +# {{{ vertex matching + +def _match_vertices( + 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: + # 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] + + mapped_src_vertices = aff_map(src_vertices) + + 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)]]) + + 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) + + 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 + +# }}} + + +# {{{ 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 _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, 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 + :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. + :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 + 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 + + 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=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]] + nunmatched = len(unmatched_bdry_m_vertex_indices) + if nunmatched > 0: + vertices = mesh.vertices[:, unmatched_bdry_m_vertex_indices] + 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" + + "\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]) + + 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 + + # 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( + 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, :]], + elements=face_ids.elements[face_index_pairs[1, :]], + faces=face_ids.faces[face_index_pairs[1, :]])) + +# }}} + + +# {{{ boundary gluing + +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 + 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). Currently + operates only on facial adjacency; any existing nodal adjacency in *mesh* is + ignored/invalidated. + + :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. 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 + 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, use_tree=use_tree) + for mapping, tol in bdry_pair_mappings_and_tols] + + 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, (mapping, _) in enumerate(bdry_pair_mappings_and_tols): + 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)): + 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, + elements=elements, + element_faces=element_faces, + 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) + + return mesh.copy( + nodal_adjacency=False, + facial_adjacency_groups=facial_adjacency_groups) + +# }}} + + # {{{ 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_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") @@ -941,8 +1293,72 @@ 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 - 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 162eb1e7c..e1320d96c 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) @@ -205,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 @@ -778,6 +891,138 @@ def test_cube_icosphere(actx_factory, order, visualize=True): # }}} +# {{{ mesh boundary gluing + +@pytest.mark.parametrize("use_tree", [False, True]) +def test_glued_mesh(use_tree): + 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, BoundaryPairMapping) + 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] + 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])) + + 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 = [ + 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 == map_lower_to_upper + assert upper_grp.aff_map == map_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]) + + 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 = [ + 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 == map_lower_to_upper + assert upper_grp.aff_map == map_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..3a80c8512 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -334,13 +334,15 @@ 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", 2, [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() @@ -396,6 +398,22 @@ def f(x): nelements_side=mesh_par, group_cls=group_cls) h = 1/mesh_par + elif mesh_name == "periodic": + assert dim == 2 or dim == 3 + + 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 + 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")