diff --git a/meshmode/discretization/connection/direct.py b/meshmode/discretization/connection/direct.py index 80d4e257b..a2672895c 100644 --- a/meshmode/discretization/connection/direct.py +++ b/meshmode/discretization/connection/direct.py @@ -250,28 +250,36 @@ def _resample_point_pick_indices(self, actx: ArrayContext, :class:`pyopencl.array.Array` containing the index subset. """ - mat = actx.to_numpy(actx.thaw( - self._resample_matrix(actx, to_group_index, ibatch_index))) - - nrows, ncols = mat.shape - result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype) + ibatch = self.groups[to_group_index].batches[ibatch_index] + from_grp = self.from_discr.groups[ibatch.from_group_index] if tol_multiplier is None: tol_multiplier = 50 - tol = np.finfo(mat.dtype).eps * tol_multiplier + tol = np.finfo(ibatch.result_unit_nodes.dtype).eps * tol_multiplier - for irow in range(nrows): - one_indices, = np.where(np.abs(mat[irow] - 1) < tol) - zero_indices, = np.where(np.abs(mat[irow]) < tol) + dim, ntgt_nodes = ibatch.result_unit_nodes.shape + if dim == 0: + assert ntgt_nodes == 1 + return actx.freeze(actx.from_numpy(np.array([0], dtype=np.int32))) - if len(one_indices) != 1: - return None - if len(zero_indices) != ncols - 1: + dist_vecs = (ibatch.result_unit_nodes.reshape(dim, -1, 1) + - from_grp.unit_nodes.reshape(dim, 1, -1)) + dists = np.sqrt(np.sum(dist_vecs**2, axis=0)) + + result = np.zeros(ntgt_nodes, dtype=self.to_discr.mesh.element_id_dtype) + + if tol_multiplier is None: + tol_multiplier = 500 + + for irow in range(ntgt_nodes): + close_indices, = np.where(dists[irow] < tol) + + if len(close_indices) != 1: return None - one_index, = one_indices - result[irow] = one_index + close_index, = close_indices + result[irow] = close_index return actx.freeze(actx.from_numpy(result)) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 4f350eddc..5a5147e96 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -193,7 +193,7 @@ def get_map_jacobian(unit_nodes): df_inv_resid[1] = 1/det * (-ata[0, 1] * atb[0] + ata[0, 0]*atb[1]) else: - # The boundary of a 3D mesh is 2D, so that's the + # The boundary of a 3D mesh is 2D, so that (the case right above) is the # highest-dimensional case we genuinely care about. # # This stinks, performance-wise, because it's not vectorized. @@ -222,6 +222,47 @@ def get_map_jacobian(unit_nodes): max_resid = np.max(np.abs(resid)) if max_resid < tol: + if 0: + near_zero = np.abs(src_unit_nodes) < 1e-5 + near_one = np.abs(src_unit_nodes-1) < 1e-5 + near_mone = np.abs(src_unit_nodes+1) < 1e-5 + reasonable = near_zero | near_one | near_mone + if np.any(~reasonable): + unreasonable_els, unreasonable_nodes = np.where( + np.any(~reasonable, axis=0)) + mapped = apply_map(src_unit_nodes) + # print(src_unit_nodes[~reasonable]) + urel = unreasonable_els[0] + tgt_urel = tgt_bdry_nodes[:, urel] + src_urel = src_bdry_nodes[:, urel] + + dist_vecs = tgt_urel.reshape(3, -1, 1) - src_urel.reshape(3, 1, -1) + dists = np.sqrt(np.sum(dist_vecs**2, axis=0)) + + mat = (dists < 1e-12) + from_indices = np.array([np.where(row)[0][0] for row in mat]) + + sunodes = src_bdry_discr.groups[i_src_grp].unit_nodes + + print("GOOD UNIT NODES") + print(sunodes[:, from_indices]) + print("BAD UNIT NODES") + print(src_unit_nodes[:, urel]) + print() + + fixed_src_unit_nodes = src_unit_nodes.copy() + fixed_src_unit_nodes[:, urel] = sunodes[:, from_indices] + mapped_fixed = apply_map(fixed_src_unit_nodes) + + print("FIXED RESIDUAL") + print(np.max(np.abs(mapped_fixed[:, urel] + - tgt_bdry_nodes[:, urel]))) + print("UNFIXED RESIDUAL") + print(np.max(np.abs(mapped[:, urel] - tgt_bdry_nodes[:, urel]))) + + + pu.db + logger.debug("_make_cross_face_batches: gauss-newton: done, " "final residual: %g", max_resid) break diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index c64d12c0e..aa3d78f75 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -1193,13 +1193,17 @@ def generate_warped_rect_mesh(dim, order, *, nelements_side=None, def m(x): result = np.empty_like(x) - result[0] = ( - 1.5*x[0] + np.cos(x[0]) - + 0.1*np.sin(10*x[1])) - result[1] = ( - 0.05*np.cos(10*x[0]) - + 1.3*x[1] + np.sin(x[1])) - if len(x) == 3: + if len(x) >= 2: + result[0] = ( + 1.5*x[0] + np.cos(x[0]) + + 0.1*np.sin(10*x[1])) + result[1] = ( + 0.05*np.cos(10*x[0]) + + 1.3*x[1] + np.sin(x[1])) + else: + result[0] = 1.5*x[0] + np.cos(x[0]) + + if len(x) >= 3: result[2] = x[2] + np.sin(x[0] / 2) / 2 return result diff --git a/test/test_connection.py b/test/test_connection.py new file mode 100644 index 000000000..baefca83e --- /dev/null +++ b/test/test_connection.py @@ -0,0 +1,109 @@ +__copyright__ = "Copyright (C) 2020 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from arraycontext import _acf # noqa: F401 +from functools import partial +import numpy as np # noqa: F401 +import numpy.linalg as la # noqa: F401 + +import meshmode # noqa: F401 +from arraycontext import ( # noqa + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests) + +from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup +from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantSimplexGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, + PolynomialRecursiveNodesGroupFactory, + ) +from meshmode.discretization import Discretization +from meshmode.discretization.connection import FACE_RESTR_ALL +import meshmode.mesh.generation as mgen + +import pytest + +import logging +logger = logging.getLogger(__name__) + + +def connection_is_permutation(actx, conn): + for i_tgrp, cgrp in enumerate(conn.groups): + for i_batch, batch in enumerate(cgrp.batches): + if not len(batch.from_element_indices): + continue + + point_pick_indices = conn._resample_point_pick_indices( + actx, i_tgrp, i_batch) + + if point_pick_indices is None: + return False + + return True + + +@pytest.mark.parametrize("group_factory", [ + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantSimplexGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, + partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), + ]) +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("order", [1, 2, 3, 4, 5]) +def test_bdry_restriction_is_permutation(actx_factory, group_factory, dim, order): + """Check that restriction to the boundary and opposite-face swap + for the element groups, orders and dimensions above is actually just + indirect access. + """ + actx = actx_factory() + + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + + mesh = mgen.generate_warped_rect_mesh(dim, order=order, n=5, + group_cls=group_cls) + + vol_discr = Discretization(actx, mesh, group_factory(order)) + from meshmode.discretization.connection import ( + make_face_restriction, make_opposite_face_connection) + bdry_connection = make_face_restriction( + actx, vol_discr, group_factory(order), + FACE_RESTR_ALL) + + assert connection_is_permutation(actx, bdry_connection) + + opp_face = make_opposite_face_connection(actx, bdry_connection) + assert connection_is_permutation(actx, opp_face) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker