From 5d6dd963a44747e1985eeb7396e1fe6aa31ff445 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 28 May 2020 12:00:24 -0500 Subject: [PATCH 1/9] add a test with a split mesh with multiple groups --- test/test_meshmode.py | 80 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 77a25cd5..d8cede25 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1272,6 +1272,86 @@ def test_is_affine_group_check(mesh_name): assert all(grp.is_affine for grp in mesh.groups) == is_affine +def split_mesh(mesh, axis, cutoff): + groups = [] + for grp in mesh.groups: + mask = np.any(mesh.vertices[axis, grp.vertex_indices] < cutoff, axis=1) + + groups.append(grp.copy( + vertex_indices=grp.vertex_indices[mask, :].copy(), + nodes=grp.nodes[:, mask, :].copy()) + ) + groups.append(grp.copy( + vertex_indices=grp.vertex_indices[~mask, :].copy(), + nodes=grp.nodes[:, ~mask, :].copy()) + ) + + from meshmode.mesh import Mesh + return Mesh( + vertices=mesh.vertices, + groups=groups, + is_conforming=mesh.is_conforming) + + +@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) + mesh = split_mesh(mesh, axis=0, cutoff=0.0) + assert mesh.facial_adjacency_groups + assert mesh.nodal_adjacency + + 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) + assert abs(cl.array.sum(bdry_f - op_bdry_f).get(queue)) < 1.0e-14 + + if boundary_tag == FACE_RESTR_ALL: + embedding = make_face_to_all_faces_embedding(conn, conn.to_discr) + check_connection(embedding) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From 949e5184b4b18b300f18c22c20e82a0216185ff2 Mon Sep 17 00:00:00 2001 From: "[6~" Date: Mon, 1 Jun 2020 15:27:33 -0500 Subject: [PATCH 2/9] Add number visualization to test_mesh_multiple_groups --- test/test_meshmode.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index d8cede25..7dab5564 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1308,6 +1308,13 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): assert mesh.facial_adjacency_groups assert mesh.nodal_adjacency + if visualize: + from meshmode.mesh.visualization import draw_2d_mesh + draw_2d_mesh(mesh, draw_element_numbers=True, draw_face_numbers=True, + set_bounding_box=True) + import matplotlib.pyplot as plt + plt.show() + from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ PolynomialWarpAndBlendGroupFactory as GroupFactory -- GitLab From f02967a14eec349d3aa2b6921c27aac17bb0aff1 Mon Sep 17 00:00:00 2001 From: "[6~" Date: Mon, 1 Jun 2020 15:29:05 -0500 Subject: [PATCH 3/9] make_opposite_face_connection: Handle multiple groups in mesh --- .../connection/opposite_face.py | 30 ++++++++++++++----- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 97570b33..1b1c9c9d 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -342,22 +342,36 @@ 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 - assert (np.array_equal( - adj.elements[adj_tgt_flags], - vbc_tgt_grp_face_batch.from_element_indices - .get(queue=queue))) + adj_els = adj.elements[adj_tgt_flags] + vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue=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(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 -- GitLab From 1a76f94fe65f3ab6b6dcb53cde645bdada244ac4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 1 Jun 2020 15:54:16 -0500 Subject: [PATCH 4/9] skip empty adjacent element batches --- .../discretization/connection/opposite_face.py | 14 ++++++++++---- test/test_meshmode.py | 7 ++++++- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 1b1c9c9d..e084525b 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -350,9 +350,14 @@ def make_opposite_face_connection(volume_to_bdry_conn): # inter-group connections. adj_tgt_flags = adj.element_faces == i_face_tgt - adj_els = adj.elements[adj_tgt_flags] - vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue=queue) + 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 @@ -362,8 +367,9 @@ def make_opposite_face_connection(volume_to_bdry_conn): 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)] + vbc_used_els = vbc_els_sort_idx[np.searchsorted( + vbc_els, adj_els, sorter=vbc_els_sort_idx + )] assert np.array_equal(vbc_els[vbc_used_els], adj_els) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 7dab5564..06b7bcfa 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1352,12 +1352,17 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): check_connection(opposite) op_bdry_f = opposite(queue, bdry_f) - assert abs(cl.array.sum(bdry_f - op_bdry_f).get(queue)) < 1.0e-14 + 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 -- GitLab From cd34cf049e3613ece24decfb17f2a0f8757bd97c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 1 Jun 2020 16:36:58 -0500 Subject: [PATCH 5/9] move split_mesh to processing and enhance a bit --- meshmode/mesh/processing.py | 45 +++++++++++++++++++++++++++++++++++++ test/test_meshmode.py | 12 ++++++++-- 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 97b9bbda..81d8526b 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -569,6 +569,51 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): # }}} +# {{{ split meshes + + +def split_mesh_groups(mesh, element_flags): + """Split all the groups in *mesh* in two according to the values of + *element_flags*. + + :arg element_flags: a :class:`numpy.ndarray` with + :attr:`~meshmode.mesh.Mesh.nelements` entries + indicating by their *Boolean* value 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*. + """ + assert element_flags.shape == (mesh.nelements,) + + new_groups = [] + for grp in mesh.groups: + mask = element_flags[ + grp.element_nr_base:grp.element_nr_base + grp.nelements] + + vertex_indices = grp.vertex_indices[mask, :].copy() + if vertex_indices.size > 0: + new_groups.append(grp.copy( + vertex_indices=vertex_indices, + nodes=grp.nodes[:, mask, :].copy() + )) + + vertex_indices = grp.vertex_indices[~mask, :].copy() + if vertex_indices.size > 0: + new_groups.append(grp.copy( + vertex_indices=vertex_indices, + nodes=grp.nodes[:, ~mask, :].copy() + )) + + from meshmode.mesh import Mesh + return Mesh( + vertices=mesh.vertices, + groups=new_groups, + is_conforming=mesh.is_conforming) + +# }}} + + # {{{ map def map_mesh(mesh, f): # noqa diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 06b7bcfa..b6a564ab 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1304,11 +1304,19 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): mesh = generate_regular_rect_mesh( a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim, n=(8,)*ambient_dim, order=order) - mesh = split_mesh(mesh, axis=0, cutoff=0.0) + 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) + mesh = split_mesh_groups(mesh, element_flags) + + assert len(mesh.groups) == 2 assert mesh.facial_adjacency_groups assert mesh.nodal_adjacency - if visualize: + if visualize and ambient_dim == 2: from meshmode.mesh.visualization import draw_2d_mesh draw_2d_mesh(mesh, draw_element_numbers=True, draw_face_numbers=True, set_bounding_box=True) -- GitLab From a4c4b0a3d6e5a951c2c09fb2a67aeb783989616c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 2 Jun 2020 14:02:19 -0500 Subject: [PATCH 6/9] make split_mesh take integer flags --- meshmode/mesh/processing.py | 54 ++++++++++++++++++++++++------------- test/test_meshmode.py | 10 ++++--- 2 files changed, 42 insertions(+), 22 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 81d8526b..e09d2d54 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 """ @@ -572,45 +573,60 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): # {{{ split meshes -def split_mesh_groups(mesh, element_flags): - """Split all the groups in *mesh* in two according to the values of - *element_flags*. +def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): + """Split all the groups in *mesh* in 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 + + .. code:: + + 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 by their *Boolean* value how the elements in a group - are to be split. + 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*. + 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 = [] - for grp in mesh.groups: - mask = element_flags[ + 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) - vertex_indices = grp.vertex_indices[mask, :].copy() - if vertex_indices.size > 0: - new_groups.append(grp.copy( - vertex_indices=vertex_indices, - nodes=grp.nodes[:, mask, :].copy() - )) + for flag in unique_grp_flags: + subgroup_to_group_map[igrp, flag] = len(new_groups) - vertex_indices = grp.vertex_indices[~mask, :].copy() - if vertex_indices.size > 0: + mask = grp_flags == flag new_groups.append(grp.copy( - vertex_indices=vertex_indices, - nodes=grp.nodes[:, ~mask, :].copy() + vertex_indices=grp.vertex_indices[mask, :].copy(), + nodes=grp.nodes[:, mask, :].copy() )) from meshmode.mesh import Mesh - return 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 + # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index b6a564ab..eb69cf7b 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1309,7 +1309,7 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): from meshmode.mesh.processing import split_mesh_groups element_flags = np.any( mesh.vertices[0, mesh.groups[0].vertex_indices] < 0.0, - axis=1) + axis=1).astype(np.int) mesh = split_mesh_groups(mesh, element_flags) assert len(mesh.groups) == 2 @@ -1318,10 +1318,14 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): if visualize and ambient_dim == 2: from meshmode.mesh.visualization import draw_2d_mesh - draw_2d_mesh(mesh, draw_element_numbers=True, draw_face_numbers=True, + 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.show() + plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300) from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ -- GitLab From 9139dcd58bddc0dd81776888462632c31a50920c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 3 Jun 2020 01:37:27 +0200 Subject: [PATCH 7/9] Apply suggestion to meshmode/mesh/processing.py --- meshmode/mesh/processing.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index e09d2d54..2ed265c0 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -577,11 +577,9 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): """Split all the groups in *mesh* in 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 + subgroups. For example, a single-group mesh with flags:: - .. code:: - - element_flags = [0, 0, 0, 42, 42, 42, 0, 0, 0, 41, 41, 41] + 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 -- GitLab From f250492a240b243f885bdade1ace9f117bcf878c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 2 Jun 2020 18:44:43 -0500 Subject: [PATCH 8/9] fix some typos --- meshmode/mesh/processing.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 2ed265c0..57300fea 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -574,7 +574,7 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): - """Split all the groups in *mesh* in according to the values of + """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:: @@ -592,7 +592,7 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): :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`` + ``(group_index, subgroup) -> new_group_index``. """ assert element_flags.shape == (mesh.nelements,) @@ -608,6 +608,7 @@ def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False): 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(), -- GitLab From 22beda11243e18d34a40c117a1665101d5863f92 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 2 Jun 2020 19:03:57 -0500 Subject: [PATCH 9/9] remove forgotten split_mesh from tests --- meshmode/mesh/processing.py | 1 - test/test_meshmode.py | 21 --------------------- 2 files changed, 22 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 57300fea..8cb8c13f 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -572,7 +572,6 @@ 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 diff --git a/test/test_meshmode.py b/test/test_meshmode.py index eb69cf7b..0031eb96 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1272,27 +1272,6 @@ def test_is_affine_group_check(mesh_name): assert all(grp.is_affine for grp in mesh.groups) == is_affine -def split_mesh(mesh, axis, cutoff): - groups = [] - for grp in mesh.groups: - mask = np.any(mesh.vertices[axis, grp.vertex_indices] < cutoff, axis=1) - - groups.append(grp.copy( - vertex_indices=grp.vertex_indices[mask, :].copy(), - nodes=grp.nodes[:, mask, :].copy()) - ) - groups.append(grp.copy( - vertex_indices=grp.vertex_indices[~mask, :].copy(), - nodes=grp.nodes[:, ~mask, :].copy()) - ) - - from meshmode.mesh import Mesh - return Mesh( - vertices=mesh.vertices, - groups=groups, - is_conforming=mesh.is_conforming) - - @pytest.mark.parametrize("ambient_dim", [1, 2, 3]) def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): ctx = ctx_factory() -- GitLab