diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 16717faee..de72491ce 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -45,6 +45,8 @@ def _make_cross_face_batches(actx, tgt_bdry_element_indices, src_bdry_element_indices, tgt_aff_map=None, src_aff_map=None): + assert len(tgt_bdry_element_indices) == len(src_bdry_element_indices) + if tgt_bdry_discr.dim == 0: return [InterpolationBatch( from_group_index=i_src_grp, @@ -431,30 +433,28 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn): # there will be separate adjacency groups for intra- and # inter-group connections. - adj_tgt_flags = adj.element_faces == i_face_tgt - adj_els = adj.elements[adj_tgt_flags] - if adj_els.size == 0: + adj_el_indices, = np.where(adj.element_faces == i_face_tgt) + if adj_el_indices.size == 0: # NOTE: this case can happen for inter-group boundaries # when all elements are adjacent on the same face # index, so all other ones will be empty continue + adj_els = adj.elements[adj_el_indices] vbc_els = thaw_to_numpy(actx, vbc_tgt_grp_face_batch.from_element_indices) if len(adj_els) == len(vbc_els): - # Same length: assert (below) that the two use the same - # ordering. vbc_used_els = slice(None) - - else: + adj_used_els = slice(None) + assert np.array_equal(vbc_els, adj_els) + elif len(vbc_els) > 0: # Genuine subset: figure out an index mapping. - vbc_els_sort_idx = np.argsort(vbc_els) - vbc_used_els = vbc_els_sort_idx[np.searchsorted( - vbc_els, adj_els, sorter=vbc_els_sort_idx - )] - - assert np.array_equal(vbc_els[vbc_used_els], adj_els) + adj_used_els, = np.where(np.in1d(adj_els, vbc_els)) + vbc_used_els, = np.where(np.in1d(vbc_els, adj_els)) + else: + # TODO: Figure out if this is actually needed + continue # find tgt_bdry_element_indices @@ -465,8 +465,10 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn): # find src_bdry_element_indices - src_vol_element_indices = adj.neighbors[adj_tgt_flags] - src_element_faces = adj.neighbor_faces[adj_tgt_flags] + src_vol_element_indices = adj.neighbors[ + adj_el_indices[adj_used_els]] + src_element_faces = adj.neighbor_faces[ + adj_el_indices[adj_used_els]] src_bdry_element_indices = src_grp_el_lookup[ src_vol_element_indices, src_element_faces] @@ -476,7 +478,8 @@ def make_opposite_face_connection(actx, volume_to_bdry_conn): # {{{ visualization (for debugging) if 0: - print("TVE", adj.elements[adj_tgt_flags]) + print("TVE", adj.elements[ + adj_el_indices[adj_used_els]]) print("TBE", tgt_bdry_element_indices) print("FVE", src_vol_element_indices) from meshmode.mesh.visualization import draw_2d_mesh diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index ca80edf83..b2a64e993 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1283,44 +1283,44 @@ def _compute_facial_adjacency_from_vertices( face_id_pairs[0].faces[is_neighbor_adj], face_id_pairs[0].elements[is_neighbor_adj]] = True - has_bdry = not np.all(face_has_neighbor) - if has_bdry: - 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] - grp_list.append( - BoundaryAdjacencyGroup( - igroup=igrp, - boundary_tag=btag, - elements=elements, - element_faces=element_faces)) - - is_untagged = ~np.any(belongs_to_bdry, axis=0) - if np.any(is_untagged): + belongs_to_bdry = np.full( + (len(boundary_tags), grp.nfaces, grp.nelements), False) + + if face_vertex_indices_to_tags is not None: + for face, element in np.ndindex(grp.nfaces, grp.nelements): + ref_fvi = grp.face_vertex_indices()[face] + fvi = frozenset(grp.vertex_indices[element, 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, face, element] = True + + for btag_idx, btag in enumerate(boundary_tags): + element_faces, elements = np.where(belongs_to_bdry[btag_idx, :, :]) + element_faces = element_faces.astype(face_id_dtype) + elements = elements.astype(element_id_dtype) + if len(elements) > 0: grp_list.append( BoundaryAdjacencyGroup( igroup=igrp, - boundary_tag=BTAG_ALL, - elements=bdry_elements, - element_faces=bdry_element_faces)) + boundary_tag=btag, + elements=elements, + element_faces=element_faces)) + + has_untagged_bdry = np.any( + ~np.any(belongs_to_bdry, axis=0) & ~face_has_neighbor) + if has_untagged_bdry: + element_faces, elements = np.where( + np.any(belongs_to_bdry, axis=0) | ~face_has_neighbor) + element_faces = element_faces.astype(face_id_dtype) + elements = elements.astype(element_id_dtype) + grp_list.append( + BoundaryAdjacencyGroup( + igroup=igrp, + boundary_tag=BTAG_ALL, + elements=elements, + element_faces=element_faces)) facial_adjacency_groups.append(grp_list)