From 11c87df5ab05db7fbdafc8dbdb0970a913bead12 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Tue, 11 Oct 2022 19:32:07 -0500 Subject: [PATCH 1/9] change tag/face mapping argument of _compute_facial_adjacency_from_vertices change the structure to something that can be constructed more efficiently --- meshmode/mesh/__init__.py | 61 ++++++++++---------- meshmode/mesh/generation.py | 107 ++++++++++++++++++------------------ meshmode/mesh/io.py | 107 +++++++++++++++++++++++++++--------- 3 files changed, 164 insertions(+), 111 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 2cbaa608c..9bd42113d 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1398,21 +1398,13 @@ def vertex_index_map_func(vertices): def _compute_facial_adjacency_from_vertices( - groups, element_id_dtype, face_id_dtype, face_vertex_indices_to_tags=None + groups, element_id_dtype, face_id_dtype, tag_to_faces=None ) -> Sequence[Sequence[FacialAdjacencyGroup]]: 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_faces is None: + tag_to_faces = [{} for grp in groups] # Match up adjacent faces according to their vertex indices @@ -1490,33 +1482,38 @@ 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_faces[igrp].items(): + # Combine known tagged faces and current boundary faces into a + # single array, lexicographically sort them, and find identical + # neighboring entries to get tags + extended_elements_and_faces = np.concatenate(( + tagged_elements_and_faces, + np.stack( + (bdry_elements, bdry_element_faces), + axis=-1))) + order = np.lexsort(extended_elements_and_faces.T) + diffs = np.diff(extended_elements_and_faces[order, :], axis=0) + match_indices, = (~np.any(diffs, axis=1)).nonzero() + # lexsort is stable, so the second entry in each match corresponds + # to the yet-to-be-tagged boundary face + face_indices = ( + order[match_indices+1] + - len(tagged_elements_and_faces)) + if len(face_indices) > 0: + elements = bdry_elements[face_indices] + element_faces = bdry_element_faces[face_indices] grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, - boundary_tag=btag, + boundary_tag=tag, elements=elements, element_faces=element_faces)) + is_tagged[face_indices] = True - is_untagged = ~np.any(belongs_to_bdry, axis=0) - if np.any(is_untagged): + if np.any(~is_tagged): grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index e25da3f05..2b26d0e23 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -1211,63 +1211,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.keys(): + # 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..eed17bd3c 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,34 @@ 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) + 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 +186,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_faces = [] for group_el_type, ngroup_elements in el_type_hist.items(): if group_el_type.dimensions != mesh_bulk_dim: @@ -204,10 +218,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 +267,49 @@ 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) + group_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] + # Combine known tagged fvis and current group fvis into a single + # array, lexicographically sort them, and find identical + # neighboring entries to get tags + if tagged_fvis.shape[1] < fvis.shape[1]: + continue + padded_fvis = np.full((fvis.shape[0], tagged_fvis.shape[1]), -1) + padded_fvis[:, :fvis.shape[1]] = fvis + extended_fvis = np.concatenate((tagged_fvis, padded_fvis)) + # Make sure vertices are in the same order + extended_fvis = np.sort(extended_fvis, axis=1) + # Lexicographically sort the face vertex indices, then diff the + # result to find tagged faces with the same vertices + order = np.lexsort(extended_fvis.T) + diffs = np.diff(extended_fvis[order, :], axis=0) + match_indices, = (~np.any(diffs, axis=1)).nonzero() + # lexsort is stable, so the second entry in each match + # corresponds to the group face + face_element_indices = ( + order[match_indices+1] + - len(tagged_fvis)) + if len(face_element_indices) > 0: + elements_and_faces = np.stack( + ( + face_element_indices, + np.full(face_element_indices.shape, fid)), + axis=-1) + if tag in group_tag_to_faces: + group_tag_to_faces[tag] = np.concatenate(( + group_tag_to_faces[tag], + elements_and_faces)) + else: + group_tag_to_faces[tag] = elements_and_faces + + tag_to_faces.append(group_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 +322,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_faces) mesh = Mesh( vertices, groups, @@ -276,7 +330,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) # }}} From c582cb77b76a3c8e51fdac67158dd824400fb460 Mon Sep 17 00:00:00 2001 From: Matt Smith Date: Wed, 12 Oct 2022 17:04:08 -0500 Subject: [PATCH 2/9] simplify numpy boolean expression Co-authored-by: Alex Fikl --- meshmode/mesh/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 9bd42113d..829aaf0ba 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1513,7 +1513,7 @@ def _compute_facial_adjacency_from_vertices( element_faces=element_faces)) is_tagged[face_indices] = True - if np.any(~is_tagged): + if not np.all(is_tagged): grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, From d9f4d67c0b1d125f39965f379695e1964290bf10 Mon Sep 17 00:00:00 2001 From: Matt Smith Date: Wed, 12 Oct 2022 17:04:53 -0500 Subject: [PATCH 3/9] omit unnecessary "keys()" Co-authored-by: Alex Fikl --- meshmode/mesh/generation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 2b26d0e23..29a28e942 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -1213,7 +1213,7 @@ def generate_box_mesh(axis_coords, order=1, *, coord_dtype=np.float64, tag_to_faces = {} - for tag in boundary_tag_to_face.keys(): + 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: From 2df711612f38949d44d5c17b2699521cdaabb3dc Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 12 Oct 2022 16:32:24 -0500 Subject: [PATCH 4/9] clean up index matching code --- meshmode/mesh/__init__.py | 50 ++++++++++++++++++++++----------------- meshmode/mesh/io.py | 24 +++++++------------ 2 files changed, 36 insertions(+), 38 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 829aaf0ba..a7bea50d7 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1345,6 +1345,29 @@ 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 of dimension ``(2, nmatches)`` containing pairs of indices into + *indices* representing entries that are the same. + """ + order = np.lexsort(indices) + diffs = np.diff(indices[:, order], axis=1) + match_indices, = (~np.any(diffs, axis=0)).nonzero() + return np.stack((order[match_indices], order[match_indices+1])) + + +def _find_matching_index_pairs(left_indices, right_indices): + """ + Return an array of dimension ``(2, nmatches)`` containing pairs of indices into + *left_indices* (row 0) and *right_indices* (row 1) representing entries that + are the same. + """ + index_pairs = _find_matching_index_pairs_merged( + np.concatenate((left_indices, right_indices), axis=1)) + index_pairs[1, :] -= left_indices.shape[1] + 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,13 +1411,8 @@ 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( @@ -1486,22 +1504,10 @@ def _compute_facial_adjacency_from_vertices( is_tagged = np.full(len(bdry_elements), False) for tag, tagged_elements_and_faces in tag_to_faces[igrp].items(): - # Combine known tagged faces and current boundary faces into a - # single array, lexicographically sort them, and find identical - # neighboring entries to get tags - extended_elements_and_faces = np.concatenate(( - tagged_elements_and_faces, - np.stack( - (bdry_elements, bdry_element_faces), - axis=-1))) - order = np.lexsort(extended_elements_and_faces.T) - diffs = np.diff(extended_elements_and_faces[order, :], axis=0) - match_indices, = (~np.any(diffs, axis=1)).nonzero() - # lexsort is stable, so the second entry in each match corresponds - # to the yet-to-be-tagged boundary face - face_indices = ( - order[match_indices+1] - - len(tagged_elements_and_faces)) + face_index_pairs = _find_matching_index_pairs( + tagged_elements_and_faces.T, + np.stack((bdry_elements, bdry_element_faces))) + face_indices = face_index_pairs[1, :] if len(face_indices) > 0: elements = bdry_elements[face_indices] element_faces = bdry_element_faces[face_indices] diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index eed17bd3c..cefd1a541 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -272,26 +272,18 @@ def get_mesh(self, return_tag_to_elements_map=False): 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] - # Combine known tagged fvis and current group fvis into a single - # array, lexicographically sort them, and find identical - # neighboring entries to get tags 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) padded_fvis[:, :fvis.shape[1]] = fvis - extended_fvis = np.concatenate((tagged_fvis, padded_fvis)) - # Make sure vertices are in the same order - extended_fvis = np.sort(extended_fvis, axis=1) - # Lexicographically sort the face vertex indices, then diff the - # result to find tagged faces with the same vertices - order = np.lexsort(extended_fvis.T) - diffs = np.diff(extended_fvis[order, :], axis=0) - match_indices, = (~np.any(diffs, axis=1)).nonzero() - # lexsort is stable, so the second entry in each match - # corresponds to the group face - face_element_indices = ( - order[match_indices+1] - - len(tagged_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: elements_and_faces = np.stack( ( From 65f7f01f37c442293b03b0b5712976fa36b7f5d0 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Wed, 12 Oct 2022 17:03:30 -0500 Subject: [PATCH 5/9] add docs for tag_to_group_faces --- meshmode/mesh/__init__.py | 18 +++++++++++++----- meshmode/mesh/io.py | 20 ++++++++++---------- 2 files changed, 23 insertions(+), 15 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index a7bea50d7..769a7fb76 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 @@ -1416,13 +1416,21 @@ def vertex_index_map_func(vertices): def _compute_facial_adjacency_from_vertices( - groups, element_id_dtype, face_id_dtype, tag_to_faces=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)`` containing + the element and face indices of each tagged face in the group. + """ if not groups: return [] - if tag_to_faces is None: - tag_to_faces = [{} for grp in groups] + if tag_to_group_faces is None: + tag_to_group_faces = [{} for grp in groups] # Match up adjacent faces according to their vertex indices @@ -1503,7 +1511,7 @@ def _compute_facial_adjacency_from_vertices( is_tagged = np.full(len(bdry_elements), False) - for tag, tagged_elements_and_faces in tag_to_faces[igrp].items(): + for tag, tagged_elements_and_faces in tag_to_group_faces[igrp].items(): face_index_pairs = _find_matching_index_pairs( tagged_elements_and_faces.T, np.stack((bdry_elements, bdry_element_faces))) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index cefd1a541..f7c003a88 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -187,7 +187,7 @@ def get_mesh(self, return_tag_to_elements_map=False): group_base_elem_nr = 0 tag_to_meshwide_elements = {} - tag_to_faces = [] + tag_to_group_faces = [] for group_el_type, ngroup_elements in el_type_hist.items(): if group_el_type.dimensions != mesh_bulk_dim: @@ -267,7 +267,7 @@ def get_mesh(self, return_tag_to_elements_map=False): group_base_elem_nr += group.nelements - group_tag_to_faces = {} + tag_to_faces = {} for tag, tagged_fvis in tag_to_fvis.items(): for fid, ref_fvi in enumerate(group.face_vertex_indices()): @@ -285,19 +285,19 @@ def get_mesh(self, return_tag_to_elements_map=False): np.sort(padded_fvis, axis=1).T) face_element_indices = face_element_index_pairs[1, :] if len(face_element_indices) > 0: - elements_and_faces = np.stack( + group_faces = np.stack( ( face_element_indices, np.full(face_element_indices.shape, fid)), axis=-1) - if tag in group_tag_to_faces: - group_tag_to_faces[tag] = np.concatenate(( - group_tag_to_faces[tag], - elements_and_faces)) + if tag in tag_to_faces: + tag_to_faces[tag] = np.concatenate(( + tag_to_faces[tag], + group_faces)) else: - group_tag_to_faces[tag] = elements_and_faces + tag_to_faces[tag] = group_faces - tag_to_faces.append(group_tag_to_faces) + tag_to_group_faces.append(tag_to_faces) for tag in tag_to_meshwide_elements.keys(): tag_to_meshwide_elements[tag] = np.array( @@ -314,7 +314,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, tag_to_faces) + groups, np.int32, np.int8, tag_to_group_faces) mesh = Mesh( vertices, groups, From 27cc6a87623dcdf9f6bd33a45b6da81a6b89b33d Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Thu, 13 Oct 2022 09:22:57 -0500 Subject: [PATCH 6/9] add some dtypes --- meshmode/mesh/io.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index f7c003a88..3a43a8294 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -172,7 +172,8 @@ def get_mesh(self, return_tag_to_elements_map=False): 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) + 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 @@ -276,7 +277,8 @@ def get_mesh(self, return_tag_to_elements_map=False): # 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) + 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( @@ -288,7 +290,9 @@ def get_mesh(self, return_tag_to_elements_map=False): group_faces = np.stack( ( face_element_indices, - np.full(face_element_indices.shape, fid)), + np.full( + face_element_indices.shape, fid, + dtype=np.int32)), axis=-1) if tag in tag_to_faces: tag_to_faces[tag] = np.concatenate(( From c97dfe743182d78b7de006b97608ac106c794cba Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Mon, 17 Oct 2022 12:51:24 -0500 Subject: [PATCH 7/9] refine documentation for index matching functions --- meshmode/mesh/__init__.py | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 769a7fb76..8656c0aa5 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1347,24 +1347,44 @@ def _concatenate_face_ids(face_ids_list): def _find_matching_index_pairs_merged(indices): """ - Return an array of dimension ``(2, nmatches)`` containing pairs of indices into - *indices* representing entries that are the same. + 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 of dimension ``(2, nmatches)`` containing pairs of indices into - *left_indices* (row 0) and *right_indices* (row 1) representing entries that - are the same. + 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. """ - index_pairs = _find_matching_index_pairs_merged( - np.concatenate((left_indices, right_indices), axis=1)) + 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 From c51a326ab092b6bc0ac369201d4e5253a1fe1974 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Mon, 17 Oct 2022 12:58:55 -0500 Subject: [PATCH 8/9] refine tag_to_group_faces documentation --- meshmode/mesh/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 8656c0aa5..95eaffcdb 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1443,8 +1443,10 @@ def _compute_facial_adjacency_from_vertices( ) -> Sequence[Sequence[FacialAdjacencyGroup]]: """ :arg tag_to_group_faces: for each group, a mapping from tag to - :class:`numpy.ndarray` of shape ``(2, nfaces)`` containing - the element and face indices of each tagged face in the group. + :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 [] From 74b53737ef8f0e561fcdda4c159aeebde65146f2 Mon Sep 17 00:00:00 2001 From: Matthew Smith Date: Mon, 17 Oct 2022 13:26:23 -0500 Subject: [PATCH 9/9] clarify naming of face index variables in boundary tag processing --- meshmode/mesh/__init__.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 95eaffcdb..b930defd8 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1534,20 +1534,20 @@ def _compute_facial_adjacency_from_vertices( is_tagged = np.full(len(bdry_elements), False) for tag, tagged_elements_and_faces in tag_to_group_faces[igrp].items(): - face_index_pairs = _find_matching_index_pairs( + tagged_and_bdry_face_index_pairs = _find_matching_index_pairs( tagged_elements_and_faces.T, np.stack((bdry_elements, bdry_element_faces))) - face_indices = face_index_pairs[1, :] - if len(face_indices) > 0: - elements = bdry_elements[face_indices] - element_faces = bdry_element_faces[face_indices] + 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=tag, elements=elements, element_faces=element_faces)) - is_tagged[face_indices] = True + is_tagged[bdry_face_indices] = True if not np.all(is_tagged): grp_list.append(