From 4bcb23148ae6e9bf14c1c5e47a9cc05e1b5bb2db Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Thu, 23 Mar 2017 14:53:05 -0500
Subject: [PATCH] Fix facial adjacency computation (reported by Ellis Hoag)
---
meshmode/mesh/__init__.py | 57 ++++++++++++++++++++++++++-----------
meshmode/mesh/processing.py | 18 +++++++++---
test/test_meshmode.py | 3 ++
3 files changed, 57 insertions(+), 21 deletions(-)
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 8c0a7962..d3c9c8cf 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -637,19 +637,21 @@ class Mesh(Record):
assert len(facial_adjacency_groups) == len(self.groups)
for fagrp_map in facial_adjacency_groups:
for fagrp in six.itervalues(fagrp_map):
- grp = self.groups[fagrp.igroup]
+ nfagrp_elements, = fagrp.elements.shape
+
+ assert fagrp.element_faces.dtype == self.face_id_dtype
+ assert fagrp.element_faces.shape == (nfagrp_elements,)
- fvi = grp.face_vertex_indices()
assert fagrp.neighbors.dtype == self.element_id_dtype
- assert fagrp.neighbors.shape == (
- grp.nelements, len(fvi))
- assert fagrp.neighbors.dtype == self.face_id_dtype
- assert fagrp.neighbor_faces.shape == (
- grp.nelements, len(fvi))
-
- is_bdry = fagrp.neighbors < 0
- assert ((1 << btag_to_index[BTAG_REALLY_ALL])
- & fagrp.neighbors[is_bdry]).all(), \
+ assert fagrp.neighbors.shape == (nfagrp_elements,)
+
+ assert fagrp.neighbor_faces.dtype == self.face_id_dtype
+ assert fagrp.neighbor_faces.shape == (nfagrp_elements,)
+
+ if fagrp.ineighbor_group is None:
+ is_bdry = fagrp.neighbors < 0
+ assert ((1 << btag_to_index[BTAG_REALLY_ALL])
+ & -fagrp.neighbors[is_bdry]).all(), \
"boundary faces without BTAG_REALLY_ALL found"
from meshmode.mesh.processing import \
@@ -868,6 +870,9 @@ def _compute_facial_adjacency_from_vertices(mesh):
face_map.setdefault(
frozenset(fvi), []).append((igrp, iel_grp, fid))
+ del igrp
+ del grp
+
# maps tuples (igrp, ineighbor_group) to number of elements
group_count = {}
for face_tuples in six.itervalues(face_map):
@@ -881,6 +886,9 @@ def _compute_facial_adjacency_from_vertices(mesh):
else:
raise RuntimeError("unexpected number of adjacent faces")
+ del face_tuples
+ del igrp
+
# {{{ build facial_adjacency_groups data structure, still empty
facial_adjacency_groups = []
@@ -895,6 +903,11 @@ def _compute_facial_adjacency_from_vertices(mesh):
neighbors = np.empty(bdry_count, dtype=mesh.element_id_dtype)
neighbor_faces = np.zeros(bdry_count, dtype=mesh.face_id_dtype)
+ # Ensure uninitialized entries get noticed
+ elements.fill(-1)
+ element_faces.fill(-1)
+ neighbor_faces.fill(-1)
+
neighbors.fill(-(
mesh.boundary_tag_bit(BTAG_ALL)
| mesh.boundary_tag_bit(BTAG_REALLY_ALL)))
@@ -915,6 +928,12 @@ def _compute_facial_adjacency_from_vertices(mesh):
neighbors = np.empty(nb_count, dtype=mesh.element_id_dtype)
neighbor_faces = np.empty(nb_count, dtype=mesh.face_id_dtype)
+ # Ensure uninitialized entries get noticed
+ elements.fill(-1)
+ element_faces.fill(-1)
+ neighbors.fill(-1)
+ neighbor_faces.fill(-1)
+
grp_map[ineighbor_group] = FacialAdjacencyGroup(
igroup=igroup,
ineighbor_group=ineighbor_group,
@@ -923,20 +942,24 @@ def _compute_facial_adjacency_from_vertices(mesh):
neighbors=neighbors,
neighbor_faces=neighbor_faces)
+ del igroup
+ del ineighbor_group
+ del grp_map
+
# }}}
# maps tuples (igrp, ineighbor_group) to number of elements filled in group
fill_count = {}
for face_tuples in six.itervalues(face_map):
if len(face_tuples) == 2:
- for (igroup, iel, iface), (inb_group, inb_el, inb_face) in [
+ for (igroup, iel, iface), (ineighbor_group, inb_el, inb_face) in [
(face_tuples[0], face_tuples[1]),
(face_tuples[1], face_tuples[0]),
]:
- idx = fill_count.get((igrp, inb_grp), 0)
- fill_count[igrp, inb_grp] = idx + 1
+ idx = fill_count.get((igroup, ineighbor_group), 0)
+ fill_count[igroup, ineighbor_group] = idx + 1
- fagrp = facial_adjacency_groups[igroup][inb_grp]
+ fagrp = facial_adjacency_groups[igroup][ineighbor_group]
fagrp.elements[idx] = iel
fagrp.element_faces[idx] = iface
fagrp.neighbors[idx] = inb_el
@@ -945,8 +968,8 @@ def _compute_facial_adjacency_from_vertices(mesh):
elif len(face_tuples) == 1:
(igroup, iel, iface), = face_tuples
- idx = fill_count.get((igrp, None), 0)
- fill_count[igrp, None] = idx + 1
+ idx = fill_count.get((igroup, None), 0)
+ fill_count[igroup, None] = idx + 1
fagrp = facial_adjacency_groups[igroup][None]
fagrp.elements[idx] = iel
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 37e4ac26..31abb9be 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -41,6 +41,8 @@ __doc__ = """
"""
+# {{{ partition_mesh
+
def partition_mesh(mesh, part_per_element, part_nr):
"""
:arg mesh: A :class:`meshmode.mesh.Mesh` to be partitioned.
@@ -141,7 +143,9 @@ def partition_mesh(mesh, part_per_element, part_nr):
from meshmode.mesh import Mesh
part_mesh = Mesh(new_vertices, new_mesh_groups)
- return (part_mesh, queried_elems)
+ return part_mesh, queried_elems
+
+# }}}
# {{{ orientations
@@ -366,15 +370,19 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
# {{{ assemble new groups list
+ nodal_adjacency = None
+ facial_adjacency_groups = None
+
if single_group:
grp_cls = None
order = None
unit_nodes = None
- nodal_adjacency = None
for mesh in meshes:
if mesh._nodal_adjacency is not None:
nodal_adjacency = False
+ if mesh._facial_adjacency_groups is not None:
+ facial_adjacency_groups = False
for group in mesh.groups:
if grp_cls is None:
@@ -400,11 +408,12 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
else:
new_groups = []
- nodal_adjacency = None
for mesh, vert_base in zip(meshes, vert_bases):
if mesh._nodal_adjacency is not None:
nodal_adjacency = False
+ if mesh._facial_adjacency_groups is not None:
+ facial_adjacency_groups = False
for group in mesh.groups:
new_vertex_indices = group.vertex_indices + vert_base
@@ -415,7 +424,8 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
from meshmode.mesh import Mesh
return Mesh(vertices, new_groups, skip_tests=skip_tests,
- nodal_adjacency=nodal_adjacency)
+ nodal_adjacency=nodal_adjacency,
+ facial_adjacency_groups=facial_adjacency_groups)
# }}}
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 84b14ec9..413a107d 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -487,6 +487,9 @@ def test_merge_and_map(ctx_getter, visualize=False):
b=np.array([5, 0, 0])[:mesh.ambient_dim])
mesh3 = merge_disjoint_meshes((mesh2, mesh))
+ mesh3.facial_adjacency_groups
+
+ mesh3.copy()
if visualize:
from meshmode.discretization import Discretization
--
GitLab