diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 97570b338f9938165004a46c59fb3ec1570e0cbe..e084525ba9a6eb5001548f432c25a0a8f861173e 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -342,22 +342,42 @@ def make_opposite_face_connection(volume_to_bdry_conn): # {{{ index wrangling - # Assert that the adjacency group and the restriction - # interpolation batch and the adjacency group have the same - # element ordering. + # The elements in the adjacency group will be a subset of + # the elements in the restriction interpolation batch: + # Imagine an inter-group boundary. The volume-to-boundary + # connection will include all faces as targets, whereas + # 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: + # 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 + + vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue) + + if len(adj_els) == len(vbc_els): + # Same length: assert (below) that the two use the same + # ordering. + vbc_used_els = slice(None) + + else: + # 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( - adj.elements[adj_tgt_flags], - vbc_tgt_grp_face_batch.from_element_indices - .get(queue=queue))) + assert np.array_equal(vbc_els[vbc_used_els], adj_els) # find to_element_indices tgt_bdry_element_indices = ( vbc_tgt_grp_face_batch.to_element_indices - .get(queue=queue)) + .get(queue=queue)[vbc_used_els]) # find from_element_indices diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 97b9bbda9c99fd0dc1f537e6cbbbaf043bc6140c..8cb8c13f3b16d4c0d3fbc9862b45ac944e226082 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -37,6 +37,7 @@ __doc__ = """ .. autofunction:: perform_flips .. autofunction:: find_bounding_box .. autofunction:: merge_disjoint_meshes +.. autofunction:: split_mesh_groups .. autofunction:: map_mesh .. autofunction:: affine_map """ @@ -569,6 +570,64 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): # }}} +# {{{ split meshes + +def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): + """Split all the groups in *mesh* according to the values of + *element_flags*. The element flags are expected to be integers + defining, for each group, how the elements are to be split into + subgroups. For example, a single-group mesh with flags:: + + element_flags = [0, 0, 0, 42, 42, 42, 0, 0, 0, 41, 41, 41] + + will create three subgroups. The integer flags need not be increasing + or contiguous and can repeat across different groups (i.e. they are + group-local). + + :arg element_flags: a :class:`numpy.ndarray` with + :attr:`~meshmode.mesh.Mesh.nelements` entries + indicating how the elements in a group are to be split. + + :returns: a :class:`~meshmode.mesh.Mesh` where each group has been split + according to flags in *element_flags*. If *return_subgroup_mapping* + is *True*, it also returns a mapping of + ``(group_index, subgroup) -> new_group_index``. + + """ + assert element_flags.shape == (mesh.nelements,) + + new_groups = [] + subgroup_to_group_map = {} + + for igrp, grp in enumerate(mesh.groups): + grp_flags = element_flags[ + grp.element_nr_base:grp.element_nr_base + grp.nelements] + unique_grp_flags = np.unique(grp_flags) + + for flag in unique_grp_flags: + subgroup_to_group_map[igrp, flag] = len(new_groups) + + # NOTE: making copies to maintain contiguity of the arrays + mask = grp_flags == flag + new_groups.append(grp.copy( + vertex_indices=grp.vertex_indices[mask, :].copy(), + nodes=grp.nodes[:, mask, :].copy() + )) + + from meshmode.mesh import Mesh + mesh = Mesh( + vertices=mesh.vertices, + groups=new_groups, + is_conforming=mesh.is_conforming) + + if return_subgroup_mapping: + return mesh, subgroup_to_group_map + else: + return mesh + +# }}} + + # {{{ map def map_mesh(mesh, f): # noqa diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 77a25cd5dd4de07e1207ca137ec5b1510d268a84..0031eb96474b887d0d39c36b99be9933661ff4d1 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1272,6 +1272,89 @@ def test_is_affine_group_check(mesh_name): assert all(grp.is_affine for grp in mesh.groups) == is_affine +@pytest.mark.parametrize("ambient_dim", [1, 2, 3]) +def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + order = 4 + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim, + n=(8,)*ambient_dim, order=order) + assert len(mesh.groups) == 1 + + from meshmode.mesh.processing import split_mesh_groups + element_flags = np.any( + mesh.vertices[0, mesh.groups[0].vertex_indices] < 0.0, + axis=1).astype(np.int) + mesh = split_mesh_groups(mesh, element_flags) + + assert len(mesh.groups) == 2 + assert mesh.facial_adjacency_groups + assert mesh.nodal_adjacency + + if visualize and ambient_dim == 2: + from meshmode.mesh.visualization import draw_2d_mesh + draw_2d_mesh(mesh, + draw_vertex_numbers=False, + draw_element_numbers=True, + draw_face_numbers=False, + set_bounding_box=True) + + import matplotlib.pyplot as plt + plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendGroupFactory as GroupFactory + discr = Discretization(ctx, mesh, GroupFactory(order)) + + if visualize: + group_id = discr.empty(queue, dtype=np.int) + for igrp, grp in enumerate(discr.groups): + group_id_view = grp.view(group_id) + group_id_view.fill(igrp) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, discr, vis_order=order) + vis.write_vtk_file("test_mesh_multiple_groups.vtu", [ + ("group_id", group_id) + ], overwrite=True) + + # check face restrictions + from meshmode.discretization.connection import ( + make_face_restriction, + make_face_to_all_faces_embedding, + make_opposite_face_connection, + check_connection) + for boundary_tag in [BTAG_ALL, FACE_RESTR_INTERIOR, FACE_RESTR_ALL]: + conn = make_face_restriction(discr, GroupFactory(order), + boundary_tag=boundary_tag, + per_face_groups=False) + check_connection(conn) + + bdry_f = conn.to_discr.empty(queue) + bdry_f.fill(1.0) + + if boundary_tag == FACE_RESTR_INTERIOR: + opposite = make_opposite_face_connection(conn) + check_connection(opposite) + + op_bdry_f = opposite(queue, bdry_f) + error = abs(cl.array.sum(bdry_f - op_bdry_f).get(queue)) + assert error < 1.0e-11, error + + if boundary_tag == FACE_RESTR_ALL: + embedding = make_face_to_all_faces_embedding(conn, conn.to_discr) + check_connection(embedding) + + em_bdry_f = embedding(queue, bdry_f) + error = abs(cl.array.sum(bdry_f - em_bdry_f).get(queue)) + assert error < 1.0e-11, error + + if __name__ == "__main__": import sys if len(sys.argv) > 1: