Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow boundaries on interior faces #285

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Allow boundaries on interior faces #285

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 19 additions & 16 deletions meshmode/discretization/connection/opposite_face.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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

Expand All @@ -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]
Expand All @@ -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
Expand Down
70 changes: 35 additions & 35 deletions meshmode/mesh/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The tagged interfaces should also be removed from the "normal" facial adjacency, to avoid accidentally computing fluxes across these boundary interfaces.

(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)

Expand Down