diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 2cbaa608c..b930defd8 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -22,7 +22,7 @@ from abc import ABC, abstractmethod from dataclasses import dataclass, replace, field -from typing import Any, ClassVar, Hashable, Optional, Tuple, Type, Sequence +from typing import Any, ClassVar, Hashable, Optional, Tuple, Type, Sequence, Mapping import numpy as np import numpy.linalg as la @@ -1345,6 +1345,49 @@ def _concatenate_face_ids(face_ids_list): faces=np.concatenate([ids.faces for ids in face_ids_list])) +def _find_matching_index_pairs_merged(indices): + """ + Return an array containing pairs of indices into *indices* representing entries + that are the same. + + Given an array *indices* of shape ``(N, nindices)`` containing integer-valued + N-tuples, returns an array of shape ``(2, nmatches)`` containing pairs of + indices into the second dimension of *indices* representing tuples that are + equal. + + Returned matches are ordered such that the second element of the pair occurs + after the first in *indices*. + """ + order = np.lexsort(indices) + diffs = np.diff(indices[:, order], axis=1) + match_indices, = (~np.any(diffs, axis=0)).nonzero() + # lexsort is stable, so the second entry in each match comes after the first + # in indices + return np.stack((order[match_indices], order[match_indices+1])) + + +def _find_matching_index_pairs(left_indices, right_indices): + """ + Return an array containing pairs of indices into *left_indices* and + *right_indices* representing entries that are the same. + + Given an array *left_indices* of shape ``(N, nindices_left)`` and an array + *right_indices* of shape ``(N, nindices_right)`` containing integer-valued + N-tuples, returns an array of shape ``(2, nmatches)`` containing pairs of + indices into the second dimension of *left_indices* (first pair element) and + *right_indices* (second pair element) representing tuples that are equal. + """ + all_indices = np.concatenate((left_indices, right_indices), axis=1) + index_pairs = _find_matching_index_pairs_merged(all_indices) + + # Indices into left_indices are the same as indices into all_indices, but need + # a conversion for right_indices + index_pairs[1, :] -= left_indices.shape[1] + assert index_pairs[1, :] >= 0 # Sanity check + + return index_pairs + + def _match_faces_by_vertices(groups, face_ids, vertex_index_map_func=None): """ Return matching faces in *face_ids* (expressed as pairs of indices into @@ -1388,31 +1431,28 @@ def vertex_index_map_func(vertices): # Normalize vertex-based "face identifiers" by sorting face_vertex_indices_increasing = np.sort(face_vertex_indices, axis=0) - # Lexicographically sort the face vertex indices, then diff the result to find - # faces with the same vertices - order = np.lexsort(face_vertex_indices_increasing) - diffs = np.diff(face_vertex_indices_increasing[:, order], axis=1) - match_indices, = (~np.any(diffs, axis=0)).nonzero() - return np.stack((order[match_indices], order[match_indices+1])) + return _find_matching_index_pairs_merged(face_vertex_indices_increasing) def _compute_facial_adjacency_from_vertices( - groups, element_id_dtype, face_id_dtype, face_vertex_indices_to_tags=None + groups: Sequence[MeshElementGroup], + element_id_dtype, + face_id_dtype, + tag_to_group_faces: Optional[Sequence[Mapping[Any, np.ndarray]]] = None ) -> Sequence[Sequence[FacialAdjacencyGroup]]: + """ + :arg tag_to_group_faces: for each group, a mapping from tag to + :class:`numpy.ndarray` of shape ``(2, nfaces)`` of tagged faces in the + group, where ``tag_to_group_faces[igrp][tag][0, :]`` contains the element + indices and ``tag_to_group_faces[igrp][tag][1, :]`` contains the reference + face indices. + """ if not groups: return [] - if face_vertex_indices_to_tags is not None: - boundary_tags = { - tag - for tags in face_vertex_indices_to_tags.values() - for tag in tags - if tags is not None} - else: - boundary_tags = set() - - boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} + if tag_to_group_faces is None: + tag_to_group_faces = [{} for grp in groups] # Match up adjacent faces according to their vertex indices @@ -1490,33 +1530,26 @@ def _compute_facial_adjacency_from_vertices( bdry_element_faces, bdry_elements = np.where(~face_has_neighbor) bdry_element_faces = bdry_element_faces.astype(face_id_dtype) bdry_elements = bdry_elements.astype(element_id_dtype) - belongs_to_bdry = np.full( - (len(boundary_tags), len(bdry_elements)), False) - - if face_vertex_indices_to_tags is not None: - for i in range(len(bdry_elements)): - ref_fvi = grp.face_vertex_indices()[bdry_element_faces[i]] - fvi = frozenset(grp.vertex_indices[bdry_elements[i], ref_fvi]) - tags = face_vertex_indices_to_tags.get(fvi, None) - if tags is not None: - for tag in tags: - btag_idx = boundary_tag_to_index[tag] - belongs_to_bdry[btag_idx, i] = True - - for btag_idx, btag in enumerate(boundary_tags): - indices, = np.where(belongs_to_bdry[btag_idx, :]) - if len(indices) > 0: - elements = bdry_elements[indices] - element_faces = bdry_element_faces[indices] + + is_tagged = np.full(len(bdry_elements), False) + + for tag, tagged_elements_and_faces in tag_to_group_faces[igrp].items(): + tagged_and_bdry_face_index_pairs = _find_matching_index_pairs( + tagged_elements_and_faces.T, + np.stack((bdry_elements, bdry_element_faces))) + bdry_face_indices = tagged_and_bdry_face_index_pairs[1, :] + if len(bdry_face_indices) > 0: + elements = bdry_elements[bdry_face_indices] + element_faces = bdry_element_faces[bdry_face_indices] grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, - boundary_tag=btag, + boundary_tag=tag, elements=elements, element_faces=element_faces)) + is_tagged[bdry_face_indices] = True - is_untagged = ~np.any(belongs_to_bdry, axis=0) - if np.any(is_untagged): + if not np.all(is_tagged): grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 4f3aa5240..98e99776d 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -1244,63 +1244,62 @@ def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64, # {{{ 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()) - nbnd_tags = len(boundary_tags) - - if nbnd_tags > 0: - vert_index_to_tuple = { - vertex_indices[itup]: itup - for itup in np.ndindex(shape)} - - for ielem in range(0, grp.nelements): - for ref_fvi in grp.face_vertex_indices(): - fvi = grp.vertex_indices[ielem, ref_fvi] - try: - fvi_tuples = [vert_index_to_tuple[i] for i in fvi] - except KeyError: - # Happens for interior faces of "X" meshes because - # midpoints aren't in vert_index_to_tuple. We don't - # care about them. - continue - - for tag in boundary_tags: - # Need to map the correct face vertices to the boundary tags - for face in boundary_tag_to_face[tag]: - if len(face) != 2: - raise ValueError( - "face identifier '%s' does not " - "consist of exactly two characters" % face) - - side, axis = face - try: - axis = axes.index(axis) - except ValueError as exc: - raise ValueError( - "unrecognized axis in face identifier " - f"'{face}'") from exc - if axis >= dim: - raise ValueError("axis in face identifier '%s' " - "does not exist in %dD" % (face, dim)) - - if side == "-": - vert_crit = 0 - elif side == "+": - vert_crit = shape[axis] - 1 - else: - raise ValueError("first character of face identifier" - " '%s' is not '+' or '-'" % face) - - if all(fvi_tuple[axis] == vert_crit - for fvi_tuple in fvi_tuples): - key = frozenset(fvi) - face_vertex_indices_to_tags.setdefault(key, - []).append(tag) + tag_to_faces = {} + + for tag in boundary_tag_to_face: + # Need to map the correct face vertices to the boundary tags + for face in boundary_tag_to_face[tag]: + if len(face) != 2: + raise ValueError("face identifier '%s' does not " + "consist of exactly two characters" % face) + + side, axis = face + try: + axis = axes.index(axis) + except ValueError as exc: + raise ValueError( + f"unrecognized axis in face identifier '{face}'") from exc + if axis >= dim: + raise ValueError("axis in face identifier '%s' does not exist in %dD" + % (face, dim)) + + if side == "-": + vert_crit = 0 + elif side == "+": + vert_crit = shape[axis] - 1 + else: + raise ValueError("first character of face identifier '%s' is not" + "'+' or '-'" % face) + + for fid, ref_fvi in enumerate(grp.face_vertex_indices()): + fvis = grp.vertex_indices[:, ref_fvi] + # Only consider faces whose vertices were all vertices of the + # original box (no midpoints, etc.) + candidate_elements, = np.where( + np.all(fvis < nvertices, axis=1)) + candidate_fvi_tuples = np.stack( + np.unravel_index( + fvis[candidate_elements, :], shape), + axis=-1) + boundary_elements = candidate_elements[ + np.all( + candidate_fvi_tuples[:, :, axis] == vert_crit, axis=1)] + boundary_faces = np.stack( + ( + boundary_elements, + np.full(boundary_elements.shape, fid)), + axis=-1) + if tag in tag_to_faces: + tag_to_faces[tag] = np.concatenate(( + tag_to_faces[tag], + boundary_faces)) + else: + tag_to_faces[tag] = boundary_faces + if tag_to_faces: from meshmode.mesh import _compute_facial_adjacency_from_vertices facial_adjacency_groups = _compute_facial_adjacency_from_vertices( - [grp], np.int32, np.int8, face_vertex_indices_to_tags) + [grp], np.int32, np.int8, [tag_to_faces]) else: facial_adjacency_groups = None diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 1f03dc0ee..3a43a8294 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -131,27 +131,12 @@ def get_mesh(self, return_tag_to_elements_map=False): # {{{ build vertex numbering - # map set of face vertex indices to list of tags associated to face - face_vertex_indices_to_tags = {} vertex_gmsh_index_to_mine = {} - for element, el_vertices in enumerate(self.element_vertices): + for el_vertices in self.element_vertices: for gmsh_vertex_nr in el_vertices: if gmsh_vertex_nr not in vertex_gmsh_index_to_mine: vertex_gmsh_index_to_mine[gmsh_vertex_nr] = \ len(vertex_gmsh_index_to_mine) - if self.tags: - el_markers = self.element_markers[element] - el_tag_indexes = ( - [self.gmsh_tag_index_to_mine[t] for t in el_markers] - if el_markers is not None else []) - # record tags of boundary dimension - el_tags = [self.tags[i][0] for i in el_tag_indexes if - self.tags[i][1] == mesh_bulk_dim - 1] - el_grp_verts = {vertex_gmsh_index_to_mine[e] for e in el_vertices} - face_vertex_indices = frozenset(el_grp_verts) - if face_vertex_indices not in face_vertex_indices_to_tags: - face_vertex_indices_to_tags[face_vertex_indices] = [] - face_vertex_indices_to_tags[face_vertex_indices] += el_tags # }}} @@ -166,6 +151,35 @@ def get_mesh(self, return_tag_to_elements_map=False): # }}} + # {{{ build map from tags to arrays of face vertex indices + + tag_to_fvis = {} + + if self.tags: + for element, el_vertices in enumerate(self.element_vertices): + el_markers = self.element_markers[element] + el_tag_indexes = ( + [self.gmsh_tag_index_to_mine[t] for t in el_markers] + if el_markers is not None else []) + # record tags of boundary dimension + el_tags = [self.tags[i][0] for i in el_tag_indexes if + self.tags[i][1] == mesh_bulk_dim - 1] + el_grp_verts = [vertex_gmsh_index_to_mine[e] for e in el_vertices] + for tag in el_tags: + fvi_list = tag_to_fvis.setdefault(tag, []) + fvi_list.append(el_grp_verts) + + for tag in tag_to_fvis.keys(): + fvi_list = tag_to_fvis[tag] + max_face_vertices = max(len(fvi) for fvi in fvi_list) + fvis = np.full( + (len(fvi_list), max_face_vertices), -1, dtype=np.int32) + for i, fvi in enumerate(fvi_list): + fvis[i, :len(fvi)] = fvi + tag_to_fvis[tag] = fvis + + # }}} + from meshmode.mesh import (Mesh, SimplexElementGroup, TensorProductElementGroup) @@ -173,7 +187,8 @@ def get_mesh(self, return_tag_to_elements_map=False): group_base_elem_nr = 0 - tag_to_elements = {} + tag_to_meshwide_elements = {} + tag_to_group_faces = [] for group_el_type, ngroup_elements in el_type_hist.items(): if group_el_type.dimensions != mesh_bulk_dim: @@ -204,10 +219,9 @@ def get_mesh(self, return_tag_to_elements_map=False): if el_markers is not None: for t in el_markers: tag = self.tags[self.gmsh_tag_index_to_mine[t]][0] - if tag not in tag_to_elements: - tag_to_elements[tag] = [group_base_elem_nr + i] - else: - tag_to_elements[tag].append(group_base_elem_nr + i) + tagged_elements = tag_to_meshwide_elements.setdefault( + tag, []) + tagged_elements.append(group_base_elem_nr + i) i += 1 @@ -254,8 +268,44 @@ def get_mesh(self, return_tag_to_elements_map=False): group_base_elem_nr += group.nelements - for tag in tag_to_elements.keys(): - tag_to_elements[tag] = np.array(tag_to_elements[tag], dtype=np.int32) + tag_to_faces = {} + + for tag, tagged_fvis in tag_to_fvis.items(): + for fid, ref_fvi in enumerate(group.face_vertex_indices()): + fvis = group.vertex_indices[:, ref_fvi] + if tagged_fvis.shape[1] < fvis.shape[1]: + # tagged_fvis does not contain any faces with the right + # number of vertices + continue + padded_fvis = np.full( + (fvis.shape[0], tagged_fvis.shape[1]), -1, dtype=np.int32) + padded_fvis[:, :fvis.shape[1]] = fvis + from meshmode.mesh import _find_matching_index_pairs + face_element_index_pairs = _find_matching_index_pairs( + # Make sure the vertices are in the same order + np.sort(tagged_fvis, axis=1).T, + np.sort(padded_fvis, axis=1).T) + face_element_indices = face_element_index_pairs[1, :] + if len(face_element_indices) > 0: + group_faces = np.stack( + ( + face_element_indices, + np.full( + face_element_indices.shape, fid, + dtype=np.int32)), + axis=-1) + if tag in tag_to_faces: + tag_to_faces[tag] = np.concatenate(( + tag_to_faces[tag], + group_faces)) + else: + tag_to_faces[tag] = group_faces + + tag_to_group_faces.append(tag_to_faces) + + for tag in tag_to_meshwide_elements.keys(): + tag_to_meshwide_elements[tag] = np.array( + tag_to_meshwide_elements[tag], dtype=np.int32) # FIXME: This is heuristic. if len(bulk_el_types) == 1: @@ -268,7 +318,7 @@ def get_mesh(self, return_tag_to_elements_map=False): if is_conforming and self.tags: from meshmode.mesh import _compute_facial_adjacency_from_vertices facial_adjacency_groups = _compute_facial_adjacency_from_vertices( - groups, np.int32, np.int8, face_vertex_indices_to_tags) + groups, np.int32, np.int8, tag_to_group_faces) mesh = Mesh( vertices, groups, @@ -276,7 +326,10 @@ def get_mesh(self, return_tag_to_elements_map=False): facial_adjacency_groups=facial_adjacency_groups, **self.mesh_construction_kwargs) - return (mesh, tag_to_elements) if return_tag_to_elements_map else mesh + return ( + (mesh, tag_to_meshwide_elements) + if return_tag_to_elements_map + else mesh) # }}}