From 0196dc7b1023e451d685e971e878bc9fba6c8c8a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 21 Apr 2020 20:08:21 -0500 Subject: [PATCH] make make_opposite_face_connection work in 1d --- .../connection/opposite_face.py | 51 +++++++++++-------- meshmode/discretization/poly_element.py | 4 +- test/test_meshmode.py | 13 ++++- 3 files changed, 43 insertions(+), 25 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 8458e959..05d55483 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -35,20 +35,34 @@ logger = logging.getLogger(__name__) # {{{ _make_cross_face_batches -def _make_cross_face_batches(queue, tgt_bdry_discr, src_bdry_discr, - i_tgt_grp, i_src_grp, - tgt_bdry_element_indices, - src_bdry_element_indices): +def _make_cross_face_batches(queue, + tgt_bdry_discr, src_bdry_discr, + i_tgt_grp, i_src_grp, + tgt_bdry_element_indices, src_bdry_element_indices): + def to_dev(ary): + return cl.array.to_device(queue, ary, array_queue=None) + + from meshmode.discretization.connection.direct import InterpolationBatch + if tgt_bdry_discr.ambient_dim == 1: + yield InterpolationBatch( + from_group_index=i_src_grp, + from_element_indices=to_dev(src_bdry_element_indices), + to_element_indices=to_dev(tgt_bdry_element_indices), + result_unit_nodes=src_bdry_discr.groups[i_src_grp].unit_nodes, + to_element_face=None) + return # FIXME: This should view-then-transfer # (but PyOpenCL doesn't do non-contiguous transfers for now). - tgt_bdry_nodes = (tgt_bdry_discr.groups[i_tgt_grp].view(tgt_bdry_discr.nodes(). - get(queue=queue))[:, tgt_bdry_element_indices]) + tgt_bdry_nodes = (tgt_bdry_discr.groups[i_tgt_grp] + .view(tgt_bdry_discr.nodes().get(queue=queue)) + [:, tgt_bdry_element_indices]) # FIXME: This should view-then-transfer # (but PyOpenCL doesn't do non-contiguous transfers for now). - src_bdry_nodes = (src_bdry_discr.groups[i_src_grp].view(src_bdry_discr.nodes(). - get(queue=queue))[:, src_bdry_element_indices]) + src_bdry_nodes = (src_bdry_discr.groups[i_src_grp] + .view(src_bdry_discr.nodes().get(queue=queue)) + [:, src_bdry_element_indices]) tol = 1e4 * np.finfo(tgt_bdry_nodes.dtype).eps @@ -222,9 +236,6 @@ def _make_cross_face_batches(queue, tgt_bdry_discr, src_bdry_discr, # {{{ find groups of src_unit_nodes - def to_dev(ary): - return cl.array.to_device(queue, ary, array_queue=None) - done_elements = np.zeros(nelements, dtype=np.bool) while True: todo_elements, = np.where(~done_elements) @@ -241,11 +252,6 @@ def _make_cross_face_batches(queue, tgt_bdry_discr, src_bdry_discr, close_els = todo_elements[unit_node_dist < tol] done_elements[close_els] = True - unit_node_dist = np.max(np.max(np.abs( - src_unit_nodes[:, todo_elements, :] - - template_unit_nodes.reshape(dim, 1, -1)), - axis=2), axis=0) - from meshmode.discretization.connection.direct import InterpolationBatch yield InterpolationBatch( from_group_index=i_src_grp, @@ -388,11 +394,12 @@ def make_opposite_face_connection(volume_to_bdry_conn): # }}} - groups[i_tgt_grp].extend(_make_cross_face_batches(queue, + batches = _make_cross_face_batches(queue, bdry_discr, bdry_discr, i_tgt_grp, i_src_grp, tgt_bdry_element_indices, - src_bdry_element_indices)) + src_bdry_element_indices) + groups[i_tgt_grp].extend(batches) from meshmode.discretization.connection import ( DirectDiscretizationConnection, DiscretizationConnectionElementGroup) @@ -485,10 +492,10 @@ def make_partition_connection(local_bdry_conn, i_local_part, local_bdry_indices = local_el_lookup[elems, faces] batches = _make_cross_face_batches(queue, - local_bdry, remote_bdry, - i_local_grp, i_remote_grp, - local_bdry_indices, - remote_bdry_indices) + local_bdry, remote_bdry, + i_local_grp, i_remote_grp, + local_bdry_indices, + remote_bdry_indices) part_batches[i_local_grp].extend(batches) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 1d25e337..7aba12a0 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -141,7 +141,9 @@ class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBa @memoize_method def _quadrature_rule(self): dims = self.mesh_el_group.dim - if dims == 1: + if dims == 0: + return mp.Quadrature(np.empty((0, 1)), np.empty((0, 1))) + elif dims == 1: return mp.LegendreGaussQuadrature(self.order) else: return mp.VioreanuRokhlinSimplexQuadrature(self.order, dims) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index c73c0292..c0be60c6 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -341,6 +341,7 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, PolynomialWarpAndBlendGroupFactory ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ + ("segment", 1, [8, 16, 32]), ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [3, 5, 7]), ("warp", 3, [3, 5]), @@ -368,7 +369,15 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, for mesh_par in mesh_pars: # {{{ get mesh - if mesh_name == "blob": + if mesh_name == "segment": + assert dim == 1 + + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh( + [np.linspace(-0.5, 0.5, mesh_par)], + order=order) + h = 1.0 / mesh_par + elif mesh_name == "blob": assert dim == 2 h = mesh_par @@ -407,7 +416,6 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, bdry_x = bdry_discr.nodes()[0].with_queue(queue) bdry_f = f(bdry_x) - bdry_f_2 = opp_face(queue, bdry_f) err = la.norm((bdry_f-bdry_f_2).get(), np.inf) @@ -1102,6 +1110,7 @@ def test_vtk_overwrite(ctx_getter): # {{{ test_mesh_to_tikz + def test_mesh_to_tikz(): from meshmode.mesh.io import generate_gmsh, FileSource -- GitLab