From ba861360af38dfdbdfd49a1f122d3b458b4e1a3c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 10:27:49 -0500 Subject: [PATCH 1/7] Rewrite DirectConnection permutation in terms of node comparison --- meshmode/discretization/connection/direct.py | 34 ++++++++++++-------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/meshmode/discretization/connection/direct.py b/meshmode/discretization/connection/direct.py index 80d4e257b..742f0a40e 100644 --- a/meshmode/discretization/connection/direct.py +++ b/meshmode/discretization/connection/direct.py @@ -22,6 +22,7 @@ import numpy as np +import numpy.linalg as la import loopy as lp from pytools import memoize_in, keyed_memoize_method @@ -250,28 +251,33 @@ 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 = la.norm(dist_vecs, axis=0, ord=2) + + result = np.zeros(ntgt_nodes, dtype=self.to_discr.mesh.element_id_dtype) + + 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)) From 8506252063b151528a5e2faaaea701b7e3c11a49 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 10:29:11 -0500 Subject: [PATCH 2/7] Generalize generate_warped_rect_mesh to >3D --- meshmode/mesh/generation.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) 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 From 7fcb686bae908b1af43629a7576b8a246b128364 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 11:04:07 -0500 Subject: [PATCH 3/7] Fix flat_norm computation in dof_array --- meshmode/dof_array.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/dof_array.py b/meshmode/dof_array.py index 1d7e60d6a..a46202b68 100644 --- a/meshmode/dof_array.py +++ b/meshmode/dof_array.py @@ -667,7 +667,7 @@ def flat_norm(ary, ord=None) -> float: actx = ary.array_context return la.norm( [ - actx.np.linalg.norm(actx.np.ravel(ary, order="A"), ord=ord) + actx.np.linalg.norm(actx.np.ravel(subary, order="A"), ord=ord) for _, subary in serialize_container(ary)], ord=ord) From 9a89743d861042173d4a4a69d1ea85ef1255520e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 12:04:13 -0500 Subject: [PATCH 4/7] Delete stray whitespace in setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 8ea25d1c6..ea9d93220 100644 --- a/setup.py +++ b/setup.py @@ -50,7 +50,7 @@ def main(): # 2019.1 is required for the Firedrake CIs, which use an very specific # version of Loopy. "loopy>=2019.1", - + "arraycontext", "recursivenodes", From edd1a6a805582df311d10b14c61380748ec415ba Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 12:04:59 -0500 Subject: [PATCH 5/7] Make separate PolynomialWarpAndBlend{2,3}DGroupFactory, deprecate PolynomialWarpAndBlendGroupFactory --- examples/moving-geometry.py | 2 +- examples/plot-connectivity.py | 4 +- examples/simple-dg.py | 5 +- experiments/refinement-playground.py | 5 +- meshmode/discretization/poly_element.py | 123 ++++++++++++++++++++++- meshmode/interop/firedrake/connection.py | 44 +++----- meshmode/mesh/visualization.py | 5 +- test/test_array.py | 6 +- test/test_mesh.py | 10 +- test/test_meshmode.py | 38 +++++-- test/test_modal.py | 13 ++- test/test_partition.py | 17 ++-- test/test_refinement.py | 17 +++- test/test_visualization.py | 9 +- 14 files changed, 220 insertions(+), 78 deletions(-) diff --git a/examples/moving-geometry.py b/examples/moving-geometry.py index 24a18635b..5af6d9fcd 100644 --- a/examples/moving-geometry.py +++ b/examples/moving-geometry.py @@ -140,7 +140,7 @@ def run(actx, *, # a bit of work when reconstructing after a time step if group_factory_name == "warp_and_blend": - group_factory_cls = poly.PolynomialWarpAndBlendGroupFactory + group_factory_cls = poly.PolynomialWarpAndBlend2DRestrictingGroupFactory unit_nodes = mp.warp_and_blend_nodes(ambient_dim - 1, mesh_order) elif group_factory_name == "quadrature": diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py index 16e815410..030b58b0d 100644 --- a/examples/plot-connectivity.py +++ b/examples/plot-connectivity.py @@ -20,10 +20,10 @@ def main(): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlend3DRestrictingGroupFactory discr = Discretization( - actx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + actx, mesh, PolynomialWarpAndBlend3DRestrictingGroupFactory(order)) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, order) diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 389328c51..8f60d9022 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -74,8 +74,9 @@ def __init__(self, actx, mesh, order): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory - self.group_factory = PolynomialWarpAndBlendGroupFactory(order=order) + PolynomialWarpAndBlend2DRestrictingGroupFactory + self.group_factory = PolynomialWarpAndBlend2DRestrictingGroupFactory( + order=order) self.volume_discr = Discretization(actx, mesh, self.group_factory) assert self.volume_discr.dim == 2 diff --git a/experiments/refinement-playground.py b/experiments/refinement-playground.py index 9d43e7188..3973ff8fa 100644 --- a/experiments/refinement-playground.py +++ b/experiments/refinement-playground.py @@ -186,9 +186,10 @@ def refine_and_generate_chart_function(mesh, filename, function): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + default_simplex_group_factory discr = Discretization( - cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + cl_ctx, mesh, default_simplex_group_factory( + base_dim=mesh.dim, order=order)) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(queue, discr, order) remove_if_exists("connectivity2.vtu") diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index ef62f5249..f4c2af3d3 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -53,7 +53,8 @@ .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup -.. autoclass:: PolynomialWarpAndBlendElementGroup +.. autoclass:: PolynomialWarpAndBlend2DRestrictingElementGroup +.. autoclass:: PolynomialWarpAndBlend3DRestrictingElementGroup .. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: PolynomialGivenNodesElementGroup @@ -80,7 +81,8 @@ .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory -.. autoclass:: PolynomialWarpAndBlendGroupFactory +.. autoclass:: PolynomialWarpAndBlend2DRestrictingGroupFactory +.. autoclass:: PolynomialWarpAndBlend3DRestrictingGroupFactory .. autoclass:: PolynomialRecursiveNodesGroupFactory .. autoclass:: PolynomialEquidistantSimplexGroupFactory .. autoclass:: PolynomialGivenNodesGroupFactory @@ -289,6 +291,16 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): Uses :func:`modepy.warp_and_blend_nodes`. """ + def __init__(self, mesh_el_group, order, index): + from warnings import warn + warn("PolynomialWarpAndBlendElementGroup is deprecated, since " + "the facial restrictions of the 3D nodes are not the 2D nodes. " + "It will go away in 2022. " + "Use PolynomialWarpAndBlend2DRestrictingElementGroup or " + "PolynomialWarpAndBlend3DRestrictingElementGroup instead.", + DeprecationWarning, stacklevel=2) + super().__init__(mesh_el_group, order, index) + @property @memoize_method def _interp_nodes(self): @@ -304,6 +316,66 @@ def _interp_nodes(self): return result +class PolynomialWarpAndBlend2DRestrictingElementGroup( + _MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. + Provides nodes in two and fewer dimensions, based on the 2D + warp-and-blend nodes and their facial restrictions. + + Uses :func:`modepy.warp_and_blend_nodes`. + """ + @property + @memoize_method + def _interp_nodes(self): + dim = self.mesh_el_group.dim + if self.order == 0: + result = mp.warp_and_blend_nodes(dim, 1) + result = np.mean(result, axis=1).reshape(-1, 1) + elif dim >= 3: + raise ValueError( + "PolynomialWarpAndBlend2DRestrictingElementGroup does not " + f"provide nodes in {dim}D") + else: + result = mp.warp_and_blend_nodes(dim, self.order) + + dim2, _ = result.shape + assert dim2 == dim + return result + + +class PolynomialWarpAndBlend3DRestrictingElementGroup( + _MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. + Provides nodes in two and fewer dimensions, based on the 3D + warp-and-blend nodes and their facial restrictions. + + Uses :func:`modepy.warp_and_blend_nodes`. + """ + @property + @memoize_method + def _interp_nodes(self): + dim = self.mesh_el_group.dim + if self.order == 0: + result = mp.warp_and_blend_nodes(dim, 1) + result = np.mean(result, axis=1).reshape(-1, 1) + else: + result = mp.warp_and_blend_nodes(3, self.order) + tol = np.finfo(result.dtype).eps * 50 + for d in range(dim, 3): + wanted = np.abs(result[d] - (-1)) < tol + result = result[:, wanted] + + result = result[:dim] + + return result + + class PolynomialRecursiveNodesElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and @@ -609,10 +681,33 @@ class QuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory): + def __init__(self, order): + from warnings import warn + warn("PolynomialWarpAndBlendGroupFactory is deprecated, since " + "the facial restrictions of the 3D nodes are not the 2D nodes. " + "It will go away in 2022. " + "Use PolynomialWarpAndBlend2DRestrictingGroupFactory or " + "PolynomialWarpAndBlend3DRestrictingGroupFactory instead.", + DeprecationWarning, stacklevel=2) + + super().__init__(order) + mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialWarpAndBlendElementGroup +class PolynomialWarpAndBlend2DRestrictingGroupFactory( + HomogeneousOrderBasedGroupFactory): + mesh_group_class = _MeshSimplexElementGroup + group_class = PolynomialWarpAndBlend2DRestrictingElementGroup + + +class PolynomialWarpAndBlend3DRestrictingGroupFactory( + HomogeneousOrderBasedGroupFactory): + mesh_group_class = _MeshSimplexElementGroup + group_class = PolynomialWarpAndBlend3DRestrictingElementGroup + + class PolynomialRecursiveNodesGroupFactory(HomogeneousOrderBasedGroupFactory): def __init__(self, order, family): super().__init__(order) @@ -674,4 +769,28 @@ class LegendreGaussLobattoTensorProductGroupFactory( # }}} +# undocumented for now, mainly for internal use +def default_simplex_group_factory(base_dim, order): + """ + :arg base_dim: The dimension of the 'base' discretization to be used. + The returned group factory will also support creating lower-dimensional + discretizations. + """ + + try: + # recursivenodes is only importable in Python 3.8 since + # it uses :func:`math.comb`, so need to check if it can + # be imported. + import recursivenodes # noqa: F401 + except ImportError: + # If it cannot be imported, use warp-and-blend nodes. + if base_dim <= 2: + return PolynomialWarpAndBlend2DRestrictingGroupFactory(order) + elif base_dim == 3: + return PolynomialWarpAndBlend3DRestrictingGroupFactory(order) + else: + raise ValueError(f"no usable set of nodes found for {base_dim}D") + + return PolynomialRecursiveNodesGroupFactory(order, family="lgl") + # vim: fdm=marker diff --git a/meshmode/interop/firedrake/connection.py b/meshmode/interop/firedrake/connection.py index fbf4516a9..1b5bc4ad1 100644 --- a/meshmode/interop/firedrake/connection.py +++ b/meshmode/interop/firedrake/connection.py @@ -39,8 +39,7 @@ from meshmode.mesh.processing import get_simplex_element_flip_matrix from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendGroupFactory, - PolynomialRecursiveNodesGroupFactory, + default_simplex_group_factory, ElementGroupFactory) from meshmode.discretization import ( Discretization, InterpolatoryElementGroupBase) @@ -623,14 +622,7 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, a :class:`~meshmode.discretization.poly_element.ElementGroupFactory` whose group class is a subclass of :class:`~meshmode.discretization.InterpolatoryElementGroupBase`. - If *None*, and :mod:`recursivenodes` can be imported, - a :class:`~meshmode.discretization.poly_element.\ -PolynomialRecursiveNodesGroupFactory` with ``"lgl"`` nodes is used. - Note that :mod:`recursivenodes` may not be importable - as it uses :func:`math.comb`, which is new in Python 3.8. - In the case that :mod:`recursivenodes` cannot be successfully - imported, a :class:`~meshmode.discretization.poly_element.\ -PolynomialWarpAndBlendGroupFactory` is used. + If *None*, and a default factory is automatically selected. :arg restrict_to_boundary: (optional) If not *None*, then must be one of the following: @@ -659,6 +651,16 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, "must be be " "'Discontinuous Lagrange', not '%s'." % ufl_elt.family()) + + # If only converting a portion of the mesh near the boundary, get + # *cells_to_use* as described in + # :func:`meshmode.interop.firedrake.mesh.import_firedrake_mesh` + cells_to_use = _get_cells_to_use(fdrake_fspace.mesh(), + restrict_to_boundary) + + mm_mesh, orient = import_firedrake_mesh(fdrake_fspace.mesh(), + cells_to_use=cells_to_use) + # Make sure grp_factory is the right type if provided, and # uses an interpolatory class. if grp_factory is not None: @@ -676,17 +678,9 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, % type(grp_factory.group_class)) # If not provided, make one else: - degree = ufl_elt.degree() - try: - # recursivenodes is only importable in Python 3.8 since - # it uses :func:`math.comb`, so need to check if it can - # be imported - import recursivenodes # noqa : F401 - family = "lgl" # L-G-Legendre - grp_factory = PolynomialRecursiveNodesGroupFactory(degree, family) - except ImportError: - # If cannot be imported, uses warp-and-blend nodes - grp_factory = PolynomialWarpAndBlendGroupFactory(degree) + grp_factory = default_simplex_group_factory( + mm_mesh.dim, ufl_elt.degree()) + # validate restrict_to_boundary, if present if restrict_to_boundary is not None: firedrake_bdy_ids = fdrake_fspace.mesh().exterior_facets.unique_markers @@ -718,15 +712,7 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None, "restrict_to_boundary. Must be an int, a tuple " "of ints, or the string 'on_boundary'.") - # If only converting a portion of the mesh near the boundary, get - # *cells_to_use* as described in - # :func:`meshmode.interop.firedrake.mesh.import_firedrake_mesh` - cells_to_use = _get_cells_to_use(fdrake_fspace.mesh(), - restrict_to_boundary) - # Create to_discr - mm_mesh, orient = import_firedrake_mesh(fdrake_fspace.mesh(), - cells_to_use=cells_to_use) to_discr = Discretization(actx, mm_mesh, grp_factory) # get firedrake unit nodes and map onto meshmode reference element diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py index 36df0a4bf..8ce555073 100644 --- a/meshmode/mesh/visualization.py +++ b/meshmode/mesh/visualization.py @@ -308,9 +308,10 @@ def vtk_visualize_mesh(actx, mesh, filename, vtk_high_order=True): order = max(mgrp.order for mgrp in mesh.groups) from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + default_simplex_group_factory from meshmode.discretization import Discretization - discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + discr = Discretization(actx, mesh, default_simplex_group_factory( + mesh.dim, order)) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, order, force_equidistant=vtk_high_order) diff --git a/test/test_array.py b/test/test_array.py index a153f62f4..91645faa9 100644 --- a/test/test_array.py +++ b/test/test_array.py @@ -35,7 +35,7 @@ with_container_arithmetic) from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import PolynomialWarpAndBlendGroupFactory +from meshmode.discretization.poly_element import default_simplex_group_factory from meshmode.dof_array import DOFArray, flat_norm from pytools.obj_array import make_obj_array @@ -67,7 +67,7 @@ def test_flatten_unflatten(actx_factory): a=(-0.5,)*ambient_dim, b=(+0.5,)*ambient_dim, n=(3,)*ambient_dim, order=1) - discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3)) + discr = Discretization(actx, mesh, default_simplex_group_factory(ambient_dim, 3)) a = np.random.randn(discr.ndofs) from meshmode.dof_array import flatten, unflatten @@ -99,7 +99,7 @@ def _get_test_containers(actx, ambient_dim=2): a=(-0.5,)*ambient_dim, b=(+0.5,)*ambient_dim, n=(3,)*ambient_dim, order=1) - discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3)) + discr = Discretization(actx, mesh, default_simplex_group_factory(ambient_dim, 3)) x = thaw(discr.nodes()[0], actx) # pylint: disable=unexpected-keyword-arg, no-value-for-parameter diff --git a/test/test_mesh.py b/test/test_mesh.py index 80c51ad86..3c1e6ac71 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -32,7 +32,7 @@ from meshmode.mesh import Mesh, SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendGroupFactory, + default_simplex_group_factory, LegendreGaussLobattoTensorProductGroupFactory, ) import meshmode.mesh.generation as mgen @@ -62,9 +62,7 @@ def test_nonequal_rect_mesh_generation(actx_factory, dim, mesh_type, order=3, mesh_type=mesh_type) from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory as GroupFactory - discr = Discretization(actx, mesh, GroupFactory(3)) + discr = Discretization(actx, mesh, default_simplex_group_factory(dim, 3)) if visualize: from meshmode.discretization.visualization import make_visualizer @@ -95,7 +93,7 @@ def test_box_mesh(actx_factory, visualize=False): actx = actx_factory() discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(7)) + default_simplex_group_factory(mesh.dim, 7)) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, 7) @@ -305,7 +303,7 @@ def test_merge_and_map(actx_factory, group_cls, visualize=False): target_unit="MM", ) - discr_grp_factory = PolynomialWarpAndBlendGroupFactory(order) + discr_grp_factory = default_simplex_group_factory(base_dim=2, order=order) else: ambient_dim = 3 mesh = mgen.generate_regular_rect_mesh( diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 7d0141f0b..e0f014923 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -34,7 +34,9 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + default_simplex_group_factory, + PolynomialWarpAndBlend2DRestrictingGroupFactory, + PolynomialWarpAndBlend3DRestrictingGroupFactory, PolynomialRecursiveNodesGroupFactory, PolynomialEquidistantSimplexGroupFactory, LegendreGaussLobattoTensorProductGroupFactory @@ -51,11 +53,24 @@ logger = logging.getLogger(__name__) +def normalize_group_factory(dim, grp_factory): + if grp_factory == "warp_and_blend": + return { + 0: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 1: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 2: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 3: PolynomialWarpAndBlend3DRestrictingGroupFactory, + }[dim] + else: + assert not isinstance(grp_factory, str) + return grp_factory + + # {{{ convergence of boundary interpolation @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + "warp_and_blend", partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), # Redundant, no information gain. @@ -86,6 +101,7 @@ def test_boundary_interpolation(actx_factory, group_factory, boundary_tag, actx = actx_factory() + group_factory = normalize_group_factory(dim, group_factory) if group_factory is LegendreGaussLobattoTensorProductGroupFactory: group_cls = TensorProductElementGroup else: @@ -184,7 +200,8 @@ def f(x): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + "warp_and_blend", + partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ @@ -201,6 +218,8 @@ def test_all_faces_interpolation(actx_factory, group_factory, actx = actx_factory() + group_factory = normalize_group_factory(dim, group_factory) + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: group_cls = TensorProductElementGroup else: @@ -306,7 +325,7 @@ def f(x): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + "warp_and_blend", LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ @@ -324,6 +343,8 @@ def test_opposite_face_interpolation(actx_factory, group_factory, logging.basicConfig(level=logging.INFO) actx = actx_factory() + group_factory = normalize_group_factory(dim, group_factory) + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: group_cls = TensorProductElementGroup else: @@ -424,7 +445,7 @@ def test_orientation_3d(actx_factory, what, mesh_gen_func, visualize=False): from meshmode.discretization import Discretization discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(3)) + default_simplex_group_factory(base_dim=3, order=3)) from pytential import bind, sym @@ -478,7 +499,7 @@ def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls, actx = actx_factory() if group_cls is SimplexElementGroup: - group_factory = PolynomialWarpAndBlendGroupFactory(mesh_order + 3) + group_factory = default_simplex_group_factory(dim, order=mesh_order + 3) elif group_cls is TensorProductElementGroup: group_factory = LegendreGaussLobattoTensorProductGroupFactory(mesh_order + 3) else: @@ -787,7 +808,8 @@ def test_mesh_multiple_groups(actx_factory, ambient_dim, visualize=False): plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300) from meshmode.discretization import Discretization - discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + grp_factory = default_simplex_group_factory(base_dim=ambient_dim, order=order) + discr = Discretization(actx, mesh, grp_factory) if visualize: group_id = discr.empty(actx, dtype=np.int32) @@ -808,7 +830,7 @@ def test_mesh_multiple_groups(actx_factory, ambient_dim, visualize=False): check_connection) for boundary_tag in [BTAG_ALL, FACE_RESTR_INTERIOR, FACE_RESTR_ALL]: conn = make_face_restriction(actx, discr, - group_factory=PolynomialWarpAndBlendGroupFactory(order), + group_factory=grp_factory, boundary_tag=boundary_tag, per_face_groups=False) check_connection(actx, conn) diff --git a/test/test_modal.py b/test/test_modal.py index 6b470b7a2..0ffca41a6 100644 --- a/test/test_modal.py +++ b/test/test_modal.py @@ -40,7 +40,8 @@ # Simplex group factories ModalSimplexGroupFactory, InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + PolynomialWarpAndBlend2DRestrictingGroupFactory, + PolynomialWarpAndBlend3DRestrictingGroupFactory, PolynomialRecursiveNodesGroupFactory, PolynomialEquidistantSimplexGroupFactory, # Tensor product group factories @@ -62,7 +63,7 @@ @pytest.mark.parametrize("nodal_group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + PolynomialWarpAndBlend2DRestrictingGroupFactory, partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), PolynomialEquidistantSimplexGroupFactory, LegendreGaussLobattoTensorProductGroupFactory, @@ -223,7 +224,7 @@ def f(x): @pytest.mark.parametrize("nodal_group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + "warp_and_blend", LegendreGaussLobattoTensorProductGroupFactory, ]) @pytest.mark.parametrize(("dim", "mesh_pars"), [ @@ -233,6 +234,12 @@ def f(x): def test_modal_truncation(actx_factory, nodal_group_factory, dim, mesh_pars): + if nodal_group_factory == "warp_and_blend": + nodal_group_factory = { + 2: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 3: PolynomialWarpAndBlend3DRestrictingGroupFactory, + }[dim] + if nodal_group_factory is LegendreGaussLobattoTensorProductGroupFactory: group_cls = TensorProductElementGroup modal_group_factory = ModalTensorProductGroupFactory diff --git a/test/test_partition.py b/test/test_partition.py index e56e423ba..0736087f7 100644 --- a/test/test_partition.py +++ b/test/test_partition.py @@ -33,8 +33,7 @@ pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests) -from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendGroupFactory) +from meshmode.discretization.poly_element import default_simplex_group_factory from meshmode.mesh import BTAG_ALL import pytest @@ -62,12 +61,12 @@ ]) def test_partition_interpolation(actx_factory, dim, mesh_pars, num_parts, num_groups, part_method): + order = 4 + np.random.seed(42) - group_factory = PolynomialWarpAndBlendGroupFactory + group_factory = default_simplex_group_factory(base_dim=dim, order=order) actx = actx_factory() - order = 4 - def f(x): return 10.*actx.np.sin(50.*x) @@ -105,7 +104,7 @@ def f(x): connected_parts.add((i_local_part, i_remote_part)) from meshmode.discretization import Discretization - vol_discrs = [Discretization(actx, part_meshes[i], group_factory(order)) + vol_discrs = [Discretization(actx, part_meshes[i], group_factory) for i in range(num_parts)] from meshmode.mesh import BTAG_PARTITION @@ -116,12 +115,12 @@ def f(x): for i_local_part, i_remote_part in connected_parts: # Mark faces within local_mesh that are connected to remote_mesh local_bdry_conn = make_face_restriction(actx, vol_discrs[i_local_part], - group_factory(order), + group_factory, BTAG_PARTITION(i_remote_part)) # Mark faces within remote_mesh that are connected to local_mesh remote_bdry_conn = make_face_restriction(actx, vol_discrs[i_remote_part], - group_factory(order), + group_factory, BTAG_PARTITION(i_local_part)) bdry_nelements = sum( @@ -342,7 +341,7 @@ def _test_mpi_boundary_swap(dim, order, num_groups): else: local_mesh = mesh_dist.receive_mesh_part() - group_factory = PolynomialWarpAndBlendGroupFactory(order) + group_factory = default_simplex_group_factory(base_dim=dim, order=order) from arraycontext import PyOpenCLArrayContext cl_ctx = cl.create_some_context() diff --git a/test/test_refinement.py b/test/test_refinement.py index 0c033da20..76ef8f059 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -39,7 +39,8 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory, + PolynomialWarpAndBlend2DRestrictingGroupFactory, + PolynomialWarpAndBlend3DRestrictingGroupFactory, PolynomialEquidistantSimplexGroupFactory, LegendreGaussLobattoTensorProductGroupFactory, GaussLegendreTensorProductGroupFactory, @@ -160,11 +161,11 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): @pytest.mark.parametrize(("refiner_cls", "group_factory"), [ (Refiner, InterpolatoryQuadratureSimplexGroupFactory), - (Refiner, PolynomialWarpAndBlendGroupFactory), + (Refiner, "warp_and_blend"), (Refiner, PolynomialEquidistantSimplexGroupFactory), (RefinerWithoutAdjacency, InterpolatoryQuadratureSimplexGroupFactory), - (RefinerWithoutAdjacency, PolynomialWarpAndBlendGroupFactory), + (RefinerWithoutAdjacency, "warp_and_blend"), (RefinerWithoutAdjacency, PolynomialEquidistantSimplexGroupFactory), (RefinerWithoutAdjacency, LegendreGaussLobattoTensorProductGroupFactory), @@ -183,10 +184,18 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): #partial(random_refine_flags, 0.4) partial(even_refine_flags, 2) ]) -# test_refinement_connection(cl._csc, RefinerWithoutAdjacency, PolynomialWarpAndBlendGroupFactory, 'warp', 2, [4, 5, 6], 5, partial(even_refine_flags, 2)) # noqa: E501 +# test_refinement_connection(cl._csc, RefinerWithoutAdjacency, PolynomialWarpAndBlend2DRestrictingGroupFactory, 'warp', 2, [4, 5, 6], 5, partial(even_refine_flags, 2)) # noqa: E501 def test_refinement_connection( actx_factory, refiner_cls, group_factory, mesh_name, dim, mesh_pars, mesh_order, refine_flags, visualize=False): + + if group_factory == "warp_and_blend": + group_factory = { + 1: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 2: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 3: PolynomialWarpAndBlend3DRestrictingGroupFactory, + }[dim] + group_cls = group_factory.mesh_group_class if issubclass(group_cls, TensorProductElementGroup): if mesh_name in ["circle", "blob"]: diff --git a/test/test_visualization.py b/test/test_visualization.py index 54470be38..da5a3a807 100644 --- a/test/test_visualization.py +++ b/test/test_visualization.py @@ -33,7 +33,7 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendGroupFactory, + default_simplex_group_factory, InterpolatoryQuadratureSimplexGroupFactory, LegendreGaussLobattoTensorProductGroupFactory, ) @@ -224,10 +224,9 @@ def test_copy_visualizer(actx_factory, ambient_dim, visualize=True): ) from meshmode.discretization import Discretization - discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(target_order)) - translated_discr = Discretization(actx, translated_mesh, - PolynomialWarpAndBlendGroupFactory(target_order)) + grp_factory = default_simplex_group_factory(ambient_dim, target_order) + discr = Discretization(actx, mesh, grp_factory) + translated_discr = Discretization(actx, translated_mesh, grp_factory) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, target_order, force_equidistant=True) From abaf862a78edd82a296071c5c1735cbf533964e6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 12:16:36 -0500 Subject: [PATCH 6/7] Add (partially passing) test that connections that should be permuations actually are --- test/test_connection.py | 128 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 test/test_connection.py diff --git a/test/test_connection.py b/test/test_connection.py new file mode 100644 index 000000000..80378595e --- /dev/null +++ b/test/test_connection.py @@ -0,0 +1,128 @@ +__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 ( + PolynomialWarpAndBlend2DRestrictingGroupFactory, + PolynomialWarpAndBlend3DRestrictingGroupFactory, + 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", [ + "warp_and_blend", + 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 == "warp_and_blend": + group_factory = { + 2: PolynomialWarpAndBlend2DRestrictingGroupFactory, + 3: PolynomialWarpAndBlend3DRestrictingGroupFactory, + }[dim] + + if group_factory is LegendreGaussLobattoTensorProductGroupFactory: + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + + mesh = mgen.generate_warped_rect_mesh(dim, order=order, nelements_side=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) + + is_lgl = group_factory is LegendreGaussLobattoTensorProductGroupFactory + + # FIXME: This should pass unconditionally + should_pass = ( + (dim == 3 and order < 2) + or (dim == 2 and not is_lgl) + or (dim == 2 and is_lgl and order < 4) + ) + + if should_pass: + opp_face = make_opposite_face_connection(actx, bdry_connection) + assert connection_is_permutation(actx, opp_face) + else: + pytest.xfail("https://github.com/inducer/meshmode/pull/105") + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker From 351b7c4e23c5db4dbd4a612ce7dfecc2cbaa48ea Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 12:33:07 -0500 Subject: [PATCH 7/7] test_mesh_multiple_groups: Test inhomogeneous polynomial degree --- meshmode/discretization/connection/opposite_face.py | 5 +++-- test/test_meshmode.py | 7 ++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 4f350eddc..154ddcd1d 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -70,7 +70,6 @@ def _make_cross_face_batches(actx, dim = src_grp.dim _, nelements, ntgt_unit_nodes = tgt_bdry_nodes.shape - assert tgt_bdry_nodes.shape == src_bdry_nodes.shape # {{{ invert face map (using Gauss-Newton) @@ -101,7 +100,9 @@ def apply_map(unit_nodes): one_deviation = np.abs(np.sum(intp_coeffs, axis=0) - 1) assert (one_deviation < tol).all(), np.max(one_deviation) - return np.einsum("fet,aef->aet", intp_coeffs, src_bdry_nodes) + mapped = np.einsum("fet,aef->aet", intp_coeffs, src_bdry_nodes) + assert tgt_bdry_nodes.shape == mapped.shape + return mapped def get_map_jacobian(unit_nodes): # unit_nodes: (dim, nelements, ntgt_unit_nodes) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index e0f014923..97de6175a 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -808,7 +808,12 @@ def test_mesh_multiple_groups(actx_factory, ambient_dim, visualize=False): plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300) from meshmode.discretization import Discretization - grp_factory = default_simplex_group_factory(base_dim=ambient_dim, order=order) + + def grp_factory(mesh_el_group, index): + return default_simplex_group_factory( + base_dim=ambient_dim, order=order + 2 if index == 0 else order + )(mesh_el_group, index) + discr = Discretization(actx, mesh, grp_factory) if visualize: