From 4b61ee2ac1bec73af634c0f71f17b4601187882a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 14 Jul 2020 16:22:20 -0500 Subject: [PATCH] Add support for Toby Isaac's recursive nodes --- doc/conf.py | 3 +- meshmode/discretization/poly_element.py | 53 +++++++++++++++++++++++++ requirements.txt | 1 + setup.py | 1 + test/test_meshmode.py | 8 +++- 5 files changed, 63 insertions(+), 3 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index 86665362..5d9f499e 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -280,5 +280,6 @@ intersphinx_mapping = { 'https://documen.tician.de/pyopencl': None, 'https://documen.tician.de/meshpy': None, 'https://documen.tician.de/modepy': None, - 'https://documen.tician.de/loopy': None + 'https://documen.tician.de/loopy': None, + 'https://tisaac.gitlab.io/recursivenodes/': None, } diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index d0f3823c..e317e311 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -40,6 +40,7 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup +.. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup @@ -49,6 +50,7 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory +.. autoclass:: PolynomialRecursiveNodesGroupFactory .. autoclass:: PolynomialEquidistantSimplexGroupFactory .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ @@ -207,6 +209,8 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): 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. + + Uses :func:`modepy.warp_and_blend_nodes`. """ @property @memoize_method @@ -223,6 +227,41 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): 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 + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. + + Supports a choice of the base *family* of 1D nodes, see the documentation + of the *family* argument to :func:`recursivenodes.recursive_nodes`. + + Requires :mod:`recursivenodes` to be installed. + + .. [Isaac20] Tobin Isaac. Recursive, parameter-free, explicitly defined + interpolation nodes for simplices. + `Arxiv preprint `__. + + .. versionadded:: 2020.2 + """ + def __init__(self, mesh_el_group, order, family, index): + super().__init__(mesh_el_group, order, index) + self.family = family + + @property + @memoize_method + def unit_nodes(self): + dim = self.mesh_el_group.dim + + from recursivenodes import recursive_nodes + result = recursive_nodes(dim, self.order, self.family, + domain="biunit").T.copy() + + dim2, nunit_nodes = result.shape + assert dim2 == dim + return result + + class PolynomialEquidistantSimplexElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and @@ -340,6 +379,20 @@ class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory): group_class = PolynomialWarpAndBlendElementGroup +class PolynomialRecursiveNodesGroupFactory(HomogeneousOrderBasedGroupFactory): + def __init__(self, order, family): + self.order = order + self.family = family + + def __call__(self, mesh_el_group, index): + if not isinstance(mesh_el_group, _MeshSimplexElementGroup): + raise TypeError("only mesh element groups of type '%s' " + "are supported" % _MeshSimplexElementGroup.__name__) + + return PolynomialRecursiveNodesElementGroup( + mesh_el_group, self.order, self.family, index) + + class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): """ .. versionadded:: 2016.1 diff --git a/requirements.txt b/requirements.txt index 4c3ff339..574e5094 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ numpy +recursivenodes git+https://gitlab.tiker.net/inducer/pytools.git#egg=pytools git+https://gitlab.tiker.net/inducer/gmsh_interop.git#egg=gmsh_interop git+https://gitlab.tiker.net/inducer/modepy.git#egg=modepy diff --git a/setup.py b/setup.py index 810d2560..1de30e57 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ def main(): "numpy", "modepy", "gmsh_interop", + "recursivenodes", "six", "pytools>=2020.3", "pytest>=2.3", diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 17773c20..3417d0aa 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from functools import partial from six.moves import range import numpy as np import numpy.linalg as la @@ -35,6 +36,7 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, + PolynomialRecursiveNodesGroupFactory, PolynomialEquidistantSimplexGroupFactory, ) from meshmode.mesh import Mesh, BTAG_ALL @@ -207,7 +209,9 @@ def test_box_boundary_tags(dim, nelem, group_factory): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlendGroupFactory, + partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), + #partial(PolynomialRecursiveNodesGroupFactory, family="gc"), ]) @pytest.mark.parametrize("boundary_tag", [ BTAG_ALL, @@ -300,7 +304,7 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag, print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-order_slack - or eoc_rec.max_error() < 1e-14) + or eoc_rec.max_error() < 3e-14) # }}} -- GitLab