From ae24076fd680ec7732205f2530a836e5a5035525 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 2 Jun 2021 10:10:52 -0500 Subject: [PATCH 01/10] add affine transformation data to facial adjacency group --- .../connection/opposite_face.py | 37 +++++++++++--- meshmode/interop/firedrake/mesh.py | 17 +++++++ meshmode/mesh/__init__.py | 50 ++++++++++++++++++- meshmode/mesh/processing.py | 49 ++++++++++++++++-- 4 files changed, 139 insertions(+), 14 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index b3f2d3a07..607ccea72 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,24 @@ 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): + if aff_transforms is not None: + mats, vecs = aff_transforms + return ( + np.einsum("edi,ien->den", mats, x) + + vecs.T.reshape(vecs.shape[1], -1, 1)) + else: + return x + + 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 +459,10 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn): vbc_tgt_grp_face_batch.to_element_indices )[vbc_used_els] + tgt_aff_transforms = ( + adj.aff_transform_mats[adj_tgt_flags, :, :], + adj.aff_transform_vecs[adj_tgt_flags, :]) + # find src_bdry_element_indices src_vol_element_indices = adj.neighbors[adj_tgt_flags] @@ -488,7 +502,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 ( @@ -560,6 +575,9 @@ def make_partition_connection(actx, *, local_bdry_conn, i_local_part, i_remote_faces = rem_ipag.element_faces[indices] i_local_vol_elems = rem_ipag.partition_neighbors[indices] i_local_faces = rem_ipag.neighbor_faces[indices] + remote_aff_transforms = ( + rem_ipag.aff_transform_mats[indices, :, :], + rem_ipag.aff_transform_vecs[indices, :]) del indices @@ -603,13 +621,18 @@ def make_partition_connection(actx, *, local_bdry_conn, i_local_part, i_remote_faces[local_indices]] assert (matched_remote_bdry_el_indices >= 0).all() + src_aff_transforms = ( + remote_aff_transforms[0][local_indices, :, :], + remote_aff_transforms[1][local_indices, :]) + # }}} 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..7b43301b3 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -334,11 +334,18 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, int_neighbors = int_neighbors[remaining_int_facs] int_neighbor_faces = int_neighbor_faces[remaining_int_facs] + # TODO: Figure out how to get ambient_dim +# from meshmode.mesh import _make_affine_identity_transforms +# int_mats, int_vecs = _make_affine_identity_transforms( +# ambient_dim, len(int_elements)) + interconnectivity_grp = FacialAdjacencyGroup(igroup=0, ineighbor_group=0, elements=int_elements, neighbors=int_neighbors, element_faces=int_element_faces, neighbor_faces=int_neighbor_faces) +# aff_transform_mats=int_mats, +# aff_transform_vecs=int_vecs) # }}} @@ -386,11 +393,18 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, new_ext_neighbor_faces)) ext_neighbors = np.concatenate((ext_neighbors, new_ext_neighbors)) + # TODO: Figure out how to get ambient_dim +# from meshmode.mesh import _make_affine_identity_transforms +# ext_mats, ext_vecs = _make_affine_identity_transforms( +# ambient_dim, len(ext_elements)) + exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor=None, elements=ext_elements, element_faces=ext_element_faces, neighbors=ext_neighbors, neighbor_faces=ext_neighbor_faces) +# aff_transform_mats=ext_mats, +# aff_transform_vecs=ext_vecs) # }}} @@ -747,6 +761,9 @@ def flip_local_face_indices(faces, elements): element_faces=new_element_faces, neighbors=fagrp.neighbors, neighbor_faces=new_neighbor_faces) +# TODO: Find out if something extra needs to be done here +# aff_transform_mats=fagrp.aff_transform_mats, +# aff_transform_vecs=fagrp.aff_transform_vecs) facial_adjacency_groups[igroup][ineighbor_group] = new_fagrp return (Mesh(vertices, [group], diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 583536863..a21d02b99 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -468,6 +468,18 @@ class FacialAdjacencyGroup(Record): Zero if ``neighbors[i]`` is negative. + .. attribute:: aff_transform_mats + + ``np.float64 [nfagrp_elements, ambient_dim, ambient_dim]``. + ``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]``. + ``aff_transform_vecs[iface, :]`` gives the vector part of the affine + mapping from face *iface* to the corresponding neighbor face. + .. automethod:: __eq__ .. automethod:: __ne__ """ @@ -481,6 +493,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,6 +568,18 @@ 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]``. + ``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]``. + ``aff_transform_vecs[iface, :]`` gives the vector part of the affine + mapping from face *iface* to the corresponding neighbor face. + .. versionadded:: 2017.1 """ @@ -1128,6 +1154,14 @@ def vertex_index_map_func(vertices): return np.stack((order[match_indices], order[match_indices+1])) +def _make_affine_identity_transforms(ambient_dim, nfaces): + mats = np.zeros((nfaces, ambient_dim, ambient_dim), dtype=np.float64) + for idim in range(ambient_dim): + mats[:, idim, idim] = 1 + vecs = np.zeros((nfaces, ambient_dim), dtype=np.float64) + return mats, vecs + + def _compute_facial_adjacency_from_vertices(groups, boundary_tags, element_id_dtype, face_id_dtype, @@ -1199,13 +1233,17 @@ def boundary_tag_bit(boundary_tag): for i_neighbor_grp in connected_groups: is_neighbor_adj = ( is_grp_adj & (face_id_pairs[1].groups == i_neighbor_grp)) + mats, vecs = _make_affine_identity_transforms( + grp.nodes.shape[0], np.count_nonzero(is_neighbor_adj)) grp_map[i_neighbor_grp] = FacialAdjacencyGroup( igroup=igrp, ineighbor_group=i_neighbor_grp, elements=face_id_pairs[0].elements[is_neighbor_adj], element_faces=face_id_pairs[0].faces[is_neighbor_adj], neighbors=face_id_pairs[1].elements[is_neighbor_adj], - neighbor_faces=face_id_pairs[1].faces[is_neighbor_adj]) + neighbor_faces=face_id_pairs[1].faces[is_neighbor_adj], + aff_transform_mats=mats, + aff_transform_vecs=vecs) face_has_neighbor[ face_id_pairs[0].faces[is_neighbor_adj], face_id_pairs[0].elements[is_neighbor_adj]] = True @@ -1229,13 +1267,17 @@ def boundary_tag_bit(boundary_tag): tag_mask |= boundary_tag_bit(tag) neighbors[i] = -((-neighbors[i]) | tag_mask) neighbor_faces = np.zeros(len(elements), dtype=face_id_dtype) + mats, vecs = _make_affine_identity_transforms( + grp.nodes.shape[0], len(elements)) grp_map[None] = FacialAdjacencyGroup( igroup=igrp, ineighbor_group=None, elements=elements, element_faces=element_faces, neighbors=neighbors, - neighbor_faces=neighbor_faces) + neighbor_faces=neighbor_faces, + aff_transform_mats=mats, + aff_transform_vecs=vecs) facial_adjacency_groups.append(grp_map) @@ -1302,6 +1344,10 @@ 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), + "aff_transform_vecs": _numpy_array_as_python( + fagrp.aff_transform_vecs), } return ",\n ".join(f"{k}={v}" for k, v in params.items()) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 695250a75..f0d440d94 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -224,6 +224,8 @@ 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, :, :] + vecs = facial_adj.aff_transform_vecs[adj_indices, :] from meshmode.mesh import FacialAdjacencyGroup local_to_local_adjacency_groups[i_part_grp][j_part_grp] =\ @@ -232,7 +234,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 +267,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,9 +340,13 @@ 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] + mats = facial_adj.aff_transform_mats[adj_indices] + vecs = facial_adj.aff_transform_vecs[adj_indices] - pairwise_adj.append(_NonLocalAdjacencyData(elements, element_faces, - neighbor_parts, global_neighbors, neighbor_faces)) + pairwise_adj.append( + _NonLocalAdjacencyData( + elements, element_faces, neighbor_parts, global_neighbors, + neighbor_faces, mats, vecs)) if pairwise_adj: nonlocal_adj_data[i_part_grp] = _NonLocalAdjacencyData( @@ -336,7 +354,9 @@ def _collect_nonlocal_adjacency_data(mesh, part_per_elem, global_elem_to_part_el 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]), + np.concatenate([adj.aff_transform_mats for adj in pairwise_adj]), + np.concatenate([adj.aff_transform_vecs for adj in pairwise_adj])) return nonlocal_adj_data @@ -452,6 +472,9 @@ 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 = np.empty((0, mesh.ambient_dim, mesh.ambient_dim), + dtype=np.float64) + vecs = np.empty((0, mesh.ambient_dim), dtype=np.float64) elif bdry is None: # Non-local adjacency only @@ -466,6 +489,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 +501,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) + from meshmode.mesh import _make_affine_identity_transforms + mats, vecs = _make_affine_identity_transforms(mesh.ambient_dim, nfaces) else: # Both; need to merge together @@ -488,6 +515,9 @@ 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 = np.empty((nfaces, mesh.ambient_dim, mesh.ambient_dim), + dtype=np.float64) + vecs = np.empty((nfaces, mesh.ambient_dim), dtype=np.float64) # Combine lists of elements/faces and sort to assist in merging combined_elements = np.concatenate((nl.elements, bdry.elements)) @@ -509,6 +539,8 @@ 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 + mats[nonlocal_indices] = nl.aff_transform_mats + vecs[nonlocal_indices] = nl.aff_transform_vecs # Merge boundary part bdry_indices = np.where(perm >= nnonlocal)[0] @@ -518,6 +550,11 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element, neighbor_parts[bdry_indices] = -1 neighbor_elements[bdry_indices] = -1 neighbor_faces[bdry_indices] = 0 + from meshmode.mesh import _make_affine_identity_transforms + bdry_mats, bdry_vecs = _make_affine_identity_transforms( + mesh.ambient_dim, nbdry) + mats[bdry_indices, :, :] = bdry_mats + vecs[bdry_indices, :] = bdry_vecs from meshmode.mesh import InterPartitionAdjacencyGroup inter_partition_adj_groups.append(InterPartitionAdjacencyGroup( @@ -525,7 +562,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 From 777922cf970e9cf9a3442f50e48b969a5ae2aeb6 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 22 Jun 2021 09:42:55 -0500 Subject: [PATCH 02/10] add mesh boundary gluing --- meshmode/mesh/processing.py | 396 ++++++++++++++++++++++++++++++++++++ 1 file changed, 396 insertions(+) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index f0d440d94..40e9f6eb2 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 @@ -1058,6 +1059,401 @@ 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, tol): + face_id_pairs_for_mapping = [] + + for btag_m, btag_n, aff_transform 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 + mapped_bdry_m_vertices = mat @ 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) + + mats_for_mapping = np.stack( + mat for _, _, (mat, _) in glued_boundary_mappings) + vecs_for_mapping = np.stack( + vec 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 + + 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 = np.empty((nfaces, mesh.ambient_dim, mesh.ambient_dim), + dtype=np.float64) + vecs = np.empty((nfaces, mesh.ambient_dim), dtype=np.float64) + + 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 + mats[indices, :, :] = old_adj.aff_transform_mats + 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 + mapping_indices = mapping_indices_grp_map[ineighbor_grp] + mats[indices, :, :] = mats_for_mapping[mapping_indices, :, :] + 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], + aff_transform_mats=old_bdry_grp.aff_transform_mats[indices, :, :], + aff_transform_vecs=old_bdry_grp.aff_transform_vecs[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, tol=1e-12): + """ + 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) 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. + :arg tol: tolerance allowed between the vertex coordinates of one boundary and + the transformed vertex coordinates of another boundary 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 in glued_boundary_mappings: + aff_transform_np = np.array(aff_transform[0]), np.array(aff_transform[1]) + glued_boundary_mappings_both_ways.append( + (btag_m, btag_n, aff_transform_np)) + if (btag_n, btag_m) not in mapped_btags: + transform_mat, transform_vec = aff_transform_np + inv_transform_mat = la.inv(transform_mat) + inv_transform_vec = -inv_transform_mat @ transform_vec + glued_boundary_mappings_both_ways.append( + (btag_n, btag_m, (inv_transform_mat, inv_transform_vec))) + + face_id_pairs_for_mapping = _match_boundary_faces(mesh, + glued_boundary_mappings_both_ways, tol) + + 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 From 3310cddad9efab930c3a754fd7653f903f428cd6 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 15 Jun 2021 16:13:03 -0500 Subject: [PATCH 03/10] make affine transformations in facial adjacency optional --- .../connection/opposite_face.py | 51 ++++-- meshmode/interop/firedrake/mesh.py | 17 -- meshmode/mesh/__init__.py | 69 ++++--- meshmode/mesh/processing.py | 170 +++++++++++++----- 4 files changed, 209 insertions(+), 98 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 607ccea72..12c263d69 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -54,13 +54,14 @@ def _make_cross_face_batches(actx, to_element_face=None)] def transform(aff_transforms, x): - if aff_transforms is not None: - mats, vecs = aff_transforms - return ( - np.einsum("edi,ien->den", mats, x) - + vecs.T.reshape(vecs.shape[1], -1, 1)) + 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: - return x + 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] @@ -459,9 +460,18 @@ 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 = ( - adj.aff_transform_mats[adj_tgt_flags, :, :], - adj.aff_transform_vecs[adj_tgt_flags, :]) + (mats, vecs) + if mats is not None or vecs is not None + else None) # find src_bdry_element_indices @@ -575,9 +585,15 @@ def make_partition_connection(actx, *, local_bdry_conn, i_local_part, i_remote_faces = rem_ipag.element_faces[indices] i_local_vol_elems = rem_ipag.partition_neighbors[indices] i_local_faces = rem_ipag.neighbor_faces[indices] - remote_aff_transforms = ( - rem_ipag.aff_transform_mats[indices, :, :], - rem_ipag.aff_transform_vecs[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 @@ -621,9 +637,18 @@ 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 = ( - remote_aff_transforms[0][local_indices, :, :], - remote_aff_transforms[1][local_indices, :]) + (mats, vecs) + if mats is not None or vecs is not None + else None) # }}} diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py index 7b43301b3..bbc9e762f 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -334,18 +334,11 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, int_neighbors = int_neighbors[remaining_int_facs] int_neighbor_faces = int_neighbor_faces[remaining_int_facs] - # TODO: Figure out how to get ambient_dim -# from meshmode.mesh import _make_affine_identity_transforms -# int_mats, int_vecs = _make_affine_identity_transforms( -# ambient_dim, len(int_elements)) - interconnectivity_grp = FacialAdjacencyGroup(igroup=0, ineighbor_group=0, elements=int_elements, neighbors=int_neighbors, element_faces=int_element_faces, neighbor_faces=int_neighbor_faces) -# aff_transform_mats=int_mats, -# aff_transform_vecs=int_vecs) # }}} @@ -393,18 +386,11 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, new_ext_neighbor_faces)) ext_neighbors = np.concatenate((ext_neighbors, new_ext_neighbors)) - # TODO: Figure out how to get ambient_dim -# from meshmode.mesh import _make_affine_identity_transforms -# ext_mats, ext_vecs = _make_affine_identity_transforms( -# ambient_dim, len(ext_elements)) - exterior_grp = FacialAdjacencyGroup(igroup=0, ineighbor=None, elements=ext_elements, element_faces=ext_element_faces, neighbors=ext_neighbors, neighbor_faces=ext_neighbor_faces) -# aff_transform_mats=ext_mats, -# aff_transform_vecs=ext_vecs) # }}} @@ -761,9 +747,6 @@ def flip_local_face_indices(faces, elements): element_faces=new_element_faces, neighbors=fagrp.neighbors, neighbor_faces=new_neighbor_faces) -# TODO: Find out if something extra needs to be done here -# aff_transform_mats=fagrp.aff_transform_mats, -# aff_transform_vecs=fagrp.aff_transform_vecs) facial_adjacency_groups[igroup][ineighbor_group] = new_fagrp return (Mesh(vertices, [group], diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index a21d02b99..4bbba72b6 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -470,13 +470,13 @@ class FacialAdjacencyGroup(Record): .. attribute:: aff_transform_mats - ``np.float64 [nfagrp_elements, ambient_dim, ambient_dim]``. + ``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]``. + ``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. @@ -484,6 +484,19 @@ class FacialAdjacencyGroup(Record): .. 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) @@ -570,19 +583,35 @@ class InterPartitionAdjacencyGroup(FacialAdjacencyGroup): .. attribute:: aff_transform_mats - ``np.float64 [nfagrp_elements, ambient_dim, ambient_dim]``. + ``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]``. + ``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) @@ -1154,14 +1183,6 @@ def vertex_index_map_func(vertices): return np.stack((order[match_indices], order[match_indices+1])) -def _make_affine_identity_transforms(ambient_dim, nfaces): - mats = np.zeros((nfaces, ambient_dim, ambient_dim), dtype=np.float64) - for idim in range(ambient_dim): - mats[:, idim, idim] = 1 - vecs = np.zeros((nfaces, ambient_dim), dtype=np.float64) - return mats, vecs - - def _compute_facial_adjacency_from_vertices(groups, boundary_tags, element_id_dtype, face_id_dtype, @@ -1233,17 +1254,13 @@ def boundary_tag_bit(boundary_tag): for i_neighbor_grp in connected_groups: is_neighbor_adj = ( is_grp_adj & (face_id_pairs[1].groups == i_neighbor_grp)) - mats, vecs = _make_affine_identity_transforms( - grp.nodes.shape[0], np.count_nonzero(is_neighbor_adj)) grp_map[i_neighbor_grp] = FacialAdjacencyGroup( igroup=igrp, ineighbor_group=i_neighbor_grp, elements=face_id_pairs[0].elements[is_neighbor_adj], element_faces=face_id_pairs[0].faces[is_neighbor_adj], neighbors=face_id_pairs[1].elements[is_neighbor_adj], - neighbor_faces=face_id_pairs[1].faces[is_neighbor_adj], - aff_transform_mats=mats, - aff_transform_vecs=vecs) + neighbor_faces=face_id_pairs[1].faces[is_neighbor_adj]) face_has_neighbor[ face_id_pairs[0].faces[is_neighbor_adj], face_id_pairs[0].elements[is_neighbor_adj]] = True @@ -1267,17 +1284,13 @@ def boundary_tag_bit(boundary_tag): tag_mask |= boundary_tag_bit(tag) neighbors[i] = -((-neighbors[i]) | tag_mask) neighbor_faces = np.zeros(len(elements), dtype=face_id_dtype) - mats, vecs = _make_affine_identity_transforms( - grp.nodes.shape[0], len(elements)) grp_map[None] = FacialAdjacencyGroup( igroup=igrp, ineighbor_group=None, elements=elements, element_faces=element_faces, neighbors=neighbors, - neighbor_faces=neighbor_faces, - aff_transform_mats=mats, - aff_transform_vecs=vecs) + neighbor_faces=neighbor_faces) facial_adjacency_groups.append(grp_map) @@ -1344,10 +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), - "aff_transform_vecs": _numpy_array_as_python( - fagrp.aff_transform_vecs), + "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/processing.py b/meshmode/mesh/processing.py index 40e9f6eb2..138ad7ab9 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -66,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): @@ -225,8 +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, :, :] - vecs = facial_adj.aff_transform_vecs[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] =\ @@ -341,8 +355,14 @@ 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] - mats = facial_adj.aff_transform_mats[adj_indices] - vecs = facial_adj.aff_transform_vecs[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) pairwise_adj.append( _NonLocalAdjacencyData( @@ -350,14 +370,36 @@ def _collect_nonlocal_adjacency_data(mesh, part_per_elem, global_elem_to_part_el 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.aff_transform_mats for adj in pairwise_adj]), - np.concatenate([adj.aff_transform_vecs for adj in pairwise_adj])) + concatenated_mats, + concatenated_vecs) return nonlocal_adj_data @@ -473,9 +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 = np.empty((0, mesh.ambient_dim, mesh.ambient_dim), - dtype=np.float64) - vecs = np.empty((0, mesh.ambient_dim), dtype=np.float64) + mats = None + vecs = None elif bdry is None: # Non-local adjacency only @@ -502,8 +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) - from meshmode.mesh import _make_affine_identity_transforms - mats, vecs = _make_affine_identity_transforms(mesh.ambient_dim, nfaces) + mats = None + vecs = None else: # Both; need to merge together @@ -516,9 +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 = np.empty((nfaces, mesh.ambient_dim, mesh.ambient_dim), - dtype=np.float64) - vecs = np.empty((nfaces, mesh.ambient_dim), dtype=np.float64) + 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)) @@ -540,8 +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 - mats[nonlocal_indices] = nl.aff_transform_mats - vecs[nonlocal_indices] = nl.aff_transform_vecs + 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] @@ -551,11 +599,6 @@ def _create_inter_partition_adjacency_groups(mesh, part_per_element, neighbor_parts[bdry_indices] = -1 neighbor_elements[bdry_indices] = -1 neighbor_faces[bdry_indices] = 0 - from meshmode.mesh import _make_affine_identity_transforms - bdry_mats, bdry_vecs = _make_affine_identity_transforms( - mesh.ambient_dim, nbdry) - mats[bdry_indices, :, :] = bdry_mats - vecs[bdry_indices, :] = bdry_vecs from meshmode.mesh import InterPartitionAdjacencyGroup inter_partition_adj_groups.append(InterPartitionAdjacencyGroup( @@ -1159,7 +1202,12 @@ def _match_boundary_faces(mesh, glued_boundary_mappings, tol): tree.insert(ivertex, bdry_n_vertex_bboxes[:, :, ivertex]) mat, vec = aff_transform - mapped_bdry_m_vertices = mat @ bdry_m_vertices + vec.reshape(-1, 1) + 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): @@ -1292,10 +1340,19 @@ def boundary_tag_bit(btag): from meshmode.mesh import _boundary_tag_bit return _boundary_tag_bit(boundary_tags, btag_to_index, btag) - mats_for_mapping = np.stack( - mat for _, _, (mat, _) in glued_boundary_mappings) - vecs_for_mapping = np.stack( - vec for _, _, (_, vec) in glued_boundary_mappings) + 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 @@ -1327,13 +1384,25 @@ def boundary_tag_bit(btag): 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 = np.empty((nfaces, mesh.ambient_dim, mesh.ambient_dim), - dtype=np.float64) - vecs = np.empty((nfaces, mesh.ambient_dim), dtype=np.float64) + + 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] @@ -1341,8 +1410,10 @@ def boundary_tag_bit(btag): element_faces[indices] = old_adj.element_faces neighbors[indices] = old_adj.neighbors neighbor_faces[indices] = old_adj.neighbor_faces - mats[indices, :, :] = old_adj.aff_transform_mats - vecs[indices, :] = old_adj.aff_transform_vecs + 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] @@ -1352,9 +1423,10 @@ def boundary_tag_bit(btag): element_faces[indices] = face_ids.faces neighbors[indices] = neighbor_face_ids.elements neighbor_faces[indices] = neighbor_face_ids.faces - mapping_indices = mapping_indices_grp_map[ineighbor_grp] - mats[indices, :, :] = mats_for_mapping[mapping_indices, :, :] - vecs[indices, :] = vecs_for_mapping[mapping_indices, :] + 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, @@ -1393,9 +1465,7 @@ def boundary_tag_bit(btag): elements=old_bdry_grp.elements[indices], element_faces=old_bdry_grp.element_faces[indices], neighbors=-flags, - neighbor_faces=old_bdry_grp.neighbor_faces[indices], - aff_transform_mats=old_bdry_grp.aff_transform_mats[indices, :, :], - aff_transform_vecs=old_bdry_grp.aff_transform_vecs[indices, :]) + neighbor_faces=old_bdry_grp.neighbor_faces[indices]) facial_adjacency_groups.append(fagrp_map) @@ -1429,13 +1499,29 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings, tol=1e-12): glued_boundary_mappings_both_ways = [] for btag_m, btag_n, aff_transform in glued_boundary_mappings: - aff_transform_np = np.array(aff_transform[0]), np.array(aff_transform[1]) + 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, aff_transform_np)) + (btag_m, btag_n, (transform_mat, transform_vec))) if (btag_n, btag_m) not in mapped_btags: - transform_mat, transform_vec = aff_transform_np - inv_transform_mat = la.inv(transform_mat) - inv_transform_vec = -inv_transform_mat @ transform_vec + 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 = ( + -transform_vec + if transform_vec is not None + else None) glued_boundary_mappings_both_ways.append( (btag_n, btag_m, (inv_transform_mat, inv_transform_vec))) From 71294afce3b3b57911710842f1e807ea52180c2e Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 17 Jun 2021 09:18:01 -0500 Subject: [PATCH 04/10] check for facial adjacency transforms in map_mesh --- meshmode/mesh/processing.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 138ad7ab9..fa58f7461 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1546,6 +1546,16 @@ 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") From e03de8a508f0f8cca65817eaf90d7eada123f5fc Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 17 Jun 2021 09:19:20 -0500 Subject: [PATCH 05/10] add facial adjacency transform handling in affine_map --- meshmode/mesh/processing.py | 69 ++++++++++++++++++++++++++++++++++--- 1 file changed, 64 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index fa58f7461..c7444aaea 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1603,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): From 1880b59f37868491c62c56765dc8a0a6402835ca Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 17 Jun 2021 13:39:16 -0500 Subject: [PATCH 06/10] set tolerance per mapping --- meshmode/mesh/processing.py | 41 ++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index c7444aaea..776ecc878 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1157,10 +1157,10 @@ def _compute_face_indices_from_mask(mask): return indices -def _match_boundary_faces(mesh, glued_boundary_mappings, tol): +def _match_boundary_faces(mesh, glued_boundary_mappings): face_id_pairs_for_mapping = [] - for btag_m, btag_n, aff_transform in glued_boundary_mappings: + 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) @@ -1327,8 +1327,8 @@ def _translate_boundary_pair_adjacency_into_group_pair_adjacency(mesh, 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)) + 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 @@ -1341,18 +1341,18 @@ def boundary_tag_bit(btag): 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]) + 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]) + 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]) + 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]) + for _, _, (_, vec), _ in glued_boundary_mappings]) from meshmode.mesh import FacialAdjacencyGroup @@ -1474,7 +1474,7 @@ def boundary_tag_bit(btag): facial_adjacency_groups=facial_adjacency_groups) -def glue_mesh_boundaries(mesh, glued_boundary_mappings, tol=1e-12): +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 @@ -1484,21 +1484,20 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings, tol=1e-12): 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) 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. - :arg tol: tolerance allowed between the vertex coordinates of one boundary and - the transformed vertex coordinates of another boundary when attempting to - match the two. + *(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) + for btag_m, btag_n, _, _ in glued_boundary_mappings) glued_boundary_mappings_both_ways = [] - for btag_m, btag_n, aff_transform in glued_boundary_mappings: + 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 @@ -1508,7 +1507,7 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings, tol=1e-12): 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))) + (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) @@ -1523,10 +1522,10 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings, tol=1e-12): if transform_vec is not None else None) glued_boundary_mappings_both_ways.append( - (btag_n, btag_m, (inv_transform_mat, inv_transform_vec))) + (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, tol) + glued_boundary_mappings_both_ways) face_id_pairs_for_group_pair, mapping_indices_for_group_pair = ( _translate_boundary_pair_adjacency_into_group_pair_adjacency( From bab998ee031b6a8f75a910a5cd873040a7c1c6e5 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 22 Jun 2021 09:36:06 -0500 Subject: [PATCH 07/10] add generate_annular_cylinder_slice_mesh and generate_periodic_annular_cylinder_slice_mesh --- meshmode/mesh/generation.py | 64 +++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) 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) From cb40a8b8fe61f9a49a9fdaf4e81a8e92fa502216 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 22 Jun 2021 09:53:22 -0500 Subject: [PATCH 08/10] add tests for facial adjacency transforms and mesh boundary gluing --- test/test_mesh.py | 255 ++++++++++++++++++++++++++++++++++++++++++ test/test_meshmode.py | 14 ++- 2 files changed, 266 insertions(+), 3 deletions(-) diff --git a/test/test_mesh.py b/test/test_mesh.py index 3c1e6ac71..02a5c27b1 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 e39b92e3d..3b0bfff7d 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") From c9fafeb5cd5a5db9dee78f2cf5ea70f4e65c5295 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 22 Jun 2021 11:27:59 -0500 Subject: [PATCH 09/10] fix bug in firedrake interop --- meshmode/interop/firedrake/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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, From f2212276dc85a14aefa09d647c721aca99d29aca Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 22 Jun 2021 19:23:01 -0500 Subject: [PATCH 10/10] Add ignore for spurious pylint warning --- meshmode/mesh/processing.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 776ecc878..c7f7da4d0 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1518,7 +1518,8 @@ def glue_mesh_boundaries(mesh, glued_boundary_mappings): else: inv_transform_mat = None inv_transform_vec = ( - -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(