diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index eb4ee95720b176255855a96460c7ce99b5fe874c..f3966e75b5523fca6a1f6c7cd33d88eef0c0302a 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -428,6 +428,9 @@ class InterPartitionAdj(): ``element_faces[i]`` is the face of ``elements[i]`` that has a neighbor. + .. attribute:: part_indices + ``part_indices[i]`` gives the partition index of the neighboring face. + .. attribute:: neighbors ``neighbors[i]`` gives the element number within the neighboring partiton @@ -449,38 +452,37 @@ class InterPartitionAdj(): self.element_faces = [] self.neighbors = [] self.neighbor_faces = [] - self.neighbor_groups = [] self.part_indices = [] - def add_connection(self, elem, face, part_idx, neighbor_group, neighbor_elem, neighbor_face): + def add_connection(self, elem, face, part_idx, neighbor_elem, neighbor_face): """ Adds a connection from ``elem`` and ``face`` within :class:`Mesh` to - ``neighbor_elem`` and ``neighbor_face`` of another neighboring partion - of type :class:`Mesh`. - :arg elem - :arg face - :arg part_idx - :arg neighbor_elem - :arg neighbor_face + ``neighbor_elem`` and ``neighbor_face`` of the neighboring partion + of type :class:`Mesh` given by `part_idx`. + :arg elem: + :arg face: + :arg part_idx: + :arg neighbor_elem: + :arg neighbor_face: """ self.elements.append(elem) self.element_faces.append(face) self.part_indices.append(part_idx) self.neighbors.append(neighbor_elem) - self.neighbor_groups.append(neighbor_group) self.neighbor_faces.append(neighbor_face) def get_neighbor(self, elem, face): """ :arg elem :arg face - :returns: A tuple ``(part_idx, neighbor_group, neighbor_elem, neighbor_face)`` of + :returns: A tuple ``(part_idx, neighbor_elem, neighbor_face)`` of neighboring elements within another :class:`Mesh`. """ for idx in range(len(self.elements)): if elem == self.elements[idx] and face == self.element_faces[idx]: - return (self.part_indices[idx], self.neighbor_groups[idx], - self.neighbors[idx], self.neighbor_faces[idx]) + return (self.part_indices[idx], + self.neighbors[idx], + self.neighbor_faces[idx]) raise RuntimeError("This face does not have a neighbor") # }}} @@ -855,6 +857,17 @@ class Mesh(Record): def __ne__(self, other): return not self.__eq__(other) + def find_igrp(self, elem): + """ + :arg elem: An element of the mesh. + :returns: The index of the group that `elem` belongs to. + """ + for igrp, grp in enumerate(self.groups): + if elem < grp.nelements: + return igrp + elem -= grp.nelements + raise RuntimeError("Could not find group with element ", elem) + # Design experience: Try not to add too many global data structures to the # mesh. Let the element groups be responsible for that at the mesh level. # diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 48effb1d1fe99a750a1d1964d47adb45b3c5a6f2..d1d57feeb460fea729b44871e6f28d5a2c9c8786 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -141,64 +141,59 @@ def partition_mesh(mesh, part_per_element, part_nr): new_nodes[group_num], unit_nodes=mesh_group.unit_nodes)) from meshmode.mesh import BTAG_ALL, BTAG_PARTITION - boundary_tags = [BTAG_PARTITION(n) for n in range(np.max(part_per_element))] + boundary_tags = [BTAG_PARTITION(n) for n in np.unique(part_per_element)] from meshmode.mesh import Mesh part_mesh = Mesh(new_vertices, new_mesh_groups, facial_adjacency_groups=None, boundary_tags=boundary_tags) - # FIXME I get errors when I try to copy part_mesh. from meshmode.mesh import InterPartitionAdj - part_mesh.interpart_adj_groups = [ - InterPartitionAdj() for _ in range(num_groups)] + interpart_grps = [InterPartitionAdj() for _ in range(len(part_mesh.groups))] - for igrp in range(num_groups): - elem_base = part_mesh.groups[igrp].element_nr_base + for igrp, grp in enumerate(part_mesh.groups): + elem_base = grp.element_nr_base boundary_adj = part_mesh.facial_adjacency_groups[igrp][None] boundary_elems = boundary_adj.elements boundary_faces = boundary_adj.element_faces - for elem_idx in range(len(boundary_elems)): - elem = boundary_elems[elem_idx] - face = boundary_faces[elem_idx] - tags = -boundary_adj.neighbors[elem_idx] + for adj_idx, elem in enumerate(boundary_elems): + face = boundary_faces[adj_idx] + tags = -boundary_adj.neighbors[adj_idx] assert tags >= 0, "Expected boundary tag in adjacency group." - parent_elem = queried_elems[elem] - parent_group_num = 0 - while parent_elem >= mesh.groups[parent_group_num].nelements: - parent_elem -= mesh.groups[parent_group_num].nelements - parent_group_num += 1 - assert parent_group_num < num_groups, "Unable to find neighbor." - parent_grp_elem_base = mesh.groups[parent_group_num].element_nr_base - parent_adj = mesh.facial_adjacency_groups[parent_group_num] - for n_grp_num, parent_facial_group in parent_adj.items(): + + parent_igrp = mesh.find_igrp(queried_elems[elem + elem_base]) + parent_elem_base = mesh.groups[parent_igrp].element_nr_base + parent_elem = queried_elems[elem + elem_base] - parent_elem_base + + parent_adj = mesh.facial_adjacency_groups[parent_igrp] + + for parent_facial_group in parent_adj.values(): for idx in np.where(parent_facial_group.elements == parent_elem)[0]: if parent_facial_group.neighbors[idx] >= 0 and \ - parent_facial_group.element_faces[idx] == face: + parent_facial_group.element_faces[idx] == face: rank_neighbor = (parent_facial_group.neighbors[idx] - + parent_grp_elem_base) + + parent_elem_base) rank_neighbor_face = parent_facial_group.neighbor_faces[idx] n_part_num = part_per_element[rank_neighbor] tags = tags & ~part_mesh.boundary_tag_bit(BTAG_ALL) tags = tags | part_mesh.boundary_tag_bit( BTAG_PARTITION(n_part_num)) - boundary_adj.neighbors[elem_idx] = -tags + boundary_adj.neighbors[adj_idx] = -tags # Find the neighbor element from the other partition n_elem = np.count_nonzero( part_per_element[:rank_neighbor] == n_part_num) - # TODO Test if this works with multiple groups - # Do I need to add the element number base? - part_mesh.interpart_adj_groups[igrp].add_connection( + interpart_grps[igrp].add_connection( elem + elem_base, face, n_part_num, - n_grp_num, n_elem, rank_neighbor_face) - return part_mesh, queried_elems + mesh = part_mesh.copy() + mesh.interpart_adj_groups = interpart_grps + return mesh, queried_elems # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cb36e4da372d5193387861697bfc6497a26d89c0..0b5f99f6f823cd2c902c1c14f9984cdc0455643e 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -95,12 +95,12 @@ def test_partition_mesh(): n = 5 num_parts = 7 from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(n, n, n)) - #TODO facial_adjacency_groups is not available from merge_disjoint_meshes. - #mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(n, n, n)) + mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(n, n, n)) + mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(n, n, n)) + mesh3 = generate_regular_rect_mesh(a=(1, 2, 2), b=(2, 3, 3), n=(n, n, n)) - #from meshmode.mesh.processing import merge_disjoint_meshes - #mesh = merge_disjoint_meshes([mesh1, mesh2]) + from meshmode.mesh.processing import merge_disjoint_meshes + mesh = merge_disjoint_meshes([mesh1, mesh2, mesh3]) adjacency_list = np.zeros((mesh.nelements,), dtype=set) for elem in range(mesh.nelements): @@ -129,33 +129,39 @@ def test_partition_mesh(): num_tags = np.zeros((num_parts,)) for part_num in range(num_parts): - (part, part_to_global) = new_meshes[part_num] + part, part_to_global = new_meshes[part_num] for grp_num, f_groups in enumerate(part.facial_adjacency_groups): f_grp = f_groups[None] for idx in range(len(f_grp.elements)): tag = -f_grp.neighbors[idx] assert tag >= 0 - elem = f_grp.elements[idx] + elem = f_grp.elements[idx] + part.groups[grp_num].element_nr_base face = f_grp.element_faces[idx] for n_part_num in range(num_parts): - (n_part, n_part_to_global) = new_meshes[n_part_num] + n_part, n_part_to_global = new_meshes[n_part_num] if tag & part.boundary_tag_bit(BTAG_PARTITION(n_part_num)) != 0: num_tags[n_part_num] += 1 - (n_part_idx, n_grp_num, n_elem, n_face) = part.\ + (i, n_elem, n_face) = part.\ interpart_adj_groups[grp_num].get_neighbor(elem, face) - assert n_part_idx == n_part_num - assert (part_num, grp_num, elem, face) == n_part.\ + assert i == n_part_num + n_grp_num = n_part.find_igrp(n_elem) + assert (part_num, elem, face) == n_part.\ interpart_adj_groups[n_grp_num].\ get_neighbor(n_elem, n_face),\ "InterpartitionAdj is not consistent" - p_elem = part_to_global[elem] + n_part_to_global = new_meshes[n_part_num][1] + p_elem = part_to_global[elem] p_n_elem = n_part_to_global[n_elem] - p_grp_num = 0 - while p_elem >= mesh.groups[p_grp_num].nelements: - p_elem -= mesh.groups[p_grp_num].nelements - p_grp_num += 1 - #p_elem_base = mesh.groups[p_grp_num].element_num_base + + p_grp_num = mesh.find_igrp(p_elem) + p_n_grp_num = mesh.find_igrp(p_n_elem) + + p_elem_base = mesh.groups[p_grp_num].element_nr_base + p_n_elem_base = mesh.groups[p_n_grp_num].element_nr_base + p_elem -= p_elem_base + p_n_elem -= p_n_elem_base + f_groups = mesh.facial_adjacency_groups[p_grp_num] for _, p_bnd_adj in f_groups.items(): for idx in range(len(p_bnd_adj.elements)):