Skip to content

Change tag/face mapping argument of _compute_facial_adjacency_from_vertices #352

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

Open
majosm wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Change tag/face mapping argument of _compute_facial_adjacency_from_vertices #352

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
11c87df
change tag/face mapping argument of _compute_facial_adjacency_from_ve…
majosm Oct 12, 2022
c582cb7
simplify numpy boolean expression
majosm Oct 12, 2022
d9f4d67
omit unnecessary "keys()"
majosm Oct 12, 2022
2df7116
clean up index matching code
majosm Oct 12, 2022
65f7f01
add docs for tag_to_group_faces
majosm Oct 12, 2022
27cc6a8
add some dtypes
majosm Oct 13, 2022
e824908
Merge remote-tracking branch 'origin/main' into optimize-tag-face-map…
majosm Oct 17, 2022
c97dfe7
refine documentation for index matching functions
majosm Oct 17, 2022
c51a326
refine tag_to_group_faces documentation
majosm Oct 17, 2022
74b5373
clarify naming of face index variables in boundary tag processing
majosm Oct 17, 2022
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
61 changes: 29 additions & 32 deletions meshmode/mesh/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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,
Expand Down
107 changes: 53 additions & 54 deletions meshmode/mesh/generation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
107 changes: 82 additions & 25 deletions meshmode/mesh/io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

# }}}

Expand All @@ -166,14 +151,43 @@ 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)

bulk_el_types = set()

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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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:
Expand All @@ -268,15 +322,18 @@ 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,
is_conforming=is_conforming,
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)

# }}}

Expand Down