diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 8c0a79628bd2268a66c5e2050f94a636388c052e..d3c9c8cfd11e4040ec7d507354776b9682f0fcae 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 37e4ac264c28f509623465c1e967fc909c070113..31abb9beb27d68090ffad264c901b01c69729b63 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 84b14ec9c735d19fe8b6e93eba36bc834465d9c2..413a107d1a3d33b52b1a0fd2eea2047822db54ed 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