diff --git a/doc/index.rst b/doc/index.rst index 7dfe6cd2..38a65a5e 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -42,6 +42,7 @@ Contents .. toctree:: :maxdepth: 2 + shapes modes nodes quadrature diff --git a/doc/modes.rst b/doc/modes.rst index 25465d69..21a75af7 100644 --- a/doc/modes.rst +++ b/doc/modes.rst @@ -1,4 +1,6 @@ Modes (Basis functions) ======================= +.. automodule:: modepy.spaces + .. automodule:: modepy.modes diff --git a/doc/nodes.rst b/doc/nodes.rst index a461abe3..c6d60997 100644 --- a/doc/nodes.rst +++ b/doc/nodes.rst @@ -1,106 +1,4 @@ Interpolation Nodes =================== -Coordinate systems on simplices -------------------------------- - -.. _tri-coords: - -Coordinates on the triangle -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Unit coordinates :math:`(r,s)`:: - - C - |\ - | \ - | O - | \ - | \ - A-----B - -Vertices in unit coordinates:: - - O = (0,0) - A = (-1,-1) - B = (1,-1) - C = (-1,1) - -Equilateral coordinates :math:`(x,y)`:: - - C - / \ - / \ - / \ - / O \ - / \ - A-----------B - -Vertices in equilateral coordinates:: - - O = (0,0) - A = (-1,-1/sqrt(3)) - B = (1,-1/sqrt(3)) - C = (0,2/sqrt(3)) - -.. _tet-coords: - -Coordinates on the tetrahedron -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Unit coordinates :math:`(r,s,t)`:: - - ^ s - | - C - /|\ - / | \ - / | \ - / | \ - / O| \ - / __A-----B---> r - /_--^ ___--^^ - ,D--^^^ - t L - -(squint, and it might start making sense...) - -Vertices in unit coordinates:: - - O=( 0, 0, 0) - A=(-1,-1,-1) - B=(+1,-1,-1) - C=(-1,+1,-1) - D=(-1,-1,+1) - -Vertices in equilateral coordinates :math:`(x,y,z)`:: - - O = (0,0,0) - A = (-1,-1/sqrt(3),-1/sqrt(6)) - B = ( 1,-1/sqrt(3),-1/sqrt(6)) - C = ( 0, 2/sqrt(3),-1/sqrt(6)) - D = ( 0, 0, 3/sqrt(6)) - -Transformations between coordinate systems ------------------------------------------- - -.. currentmodule:: modepy.tools - -All of these expect and return arrays of shape *(dims, npts)*. - -.. autofunction:: equilateral_to_unit -.. autofunction:: barycentric_to_unit -.. autofunction:: unit_to_barycentric -.. autofunction:: barycentric_to_equilateral - -Node sets for interpolation ---------------------------- - -.. currentmodule:: modepy - -.. autofunction:: equidistant_nodes -.. autofunction:: warp_and_blend_nodes -.. autofunction:: tensor_product_nodes - -Also see :class:`modepy.VioreanuRokhlinSimplexQuadrature` if nodes on the -boundary are not required. +.. automodule:: modepy.nodes diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 008bf447..dfc8eea1 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -6,10 +6,6 @@ Base classes .. automodule:: modepy.quadrature -.. currentmodule:: modepy - -.. autoclass:: Quadrature - Jacobi-Gauss quadrature in one dimension ---------------------------------------- @@ -55,4 +51,15 @@ Quadratures on the simplex .. autoclass:: VioreanuRokhlinSimplexQuadrature + +Quadratures on the hypercube +---------------------------- + +.. currentmodule:: modepy + +.. autoclass:: WitherdenVincentQuadrature + +.. autoclass:: TensorProductQuadrature +.. autoclass:: LegendreGaussTensorProductQuadrature + .. vim: sw=4 diff --git a/doc/shapes.rst b/doc/shapes.rst new file mode 100644 index 00000000..453642e5 --- /dev/null +++ b/doc/shapes.rst @@ -0,0 +1,4 @@ +Shapes +====== + +.. automodule:: modepy.shapes diff --git a/doc/tools.rst b/doc/tools.rst index b5b00069..85d1a5c1 100644 --- a/doc/tools.rst +++ b/doc/tools.rst @@ -24,11 +24,6 @@ Modal decay/residual .. automodule:: modepy.modal_decay -Interpolation quality ---------------------- - -.. currentmodule:: modepy.tools - -.. autofunction:: estimate_lebesgue_constant +.. automodule:: modepy.tools .. vim: sw=4 diff --git a/modepy/__init__.py b/modepy/__init__.py index 64228b8c..50832c47 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -22,6 +22,14 @@ """ +from modepy.shapes import ( + Shape, Face, Simplex, Hypercube, + + biunit_vertices_for_shape, faces_for_shape, + submesh_for_shape, + ) +from modepy.spaces import ( + FunctionSpace, PN, QN, space_for_shape) from modepy.modes import ( jacobi, grad_jacobi, simplex_onb, grad_simplex_onb, simplex_onb_with_mode_ids, @@ -29,19 +37,32 @@ simplex_monomial_basis_with_mode_ids, simplex_best_available_basis, grad_simplex_best_available_basis, tensor_product_basis, grad_tensor_product_basis, - legendre_tensor_product_basis, grad_legendre_tensor_product_basis) + legendre_tensor_product_basis, grad_legendre_tensor_product_basis, + symbolicize_function, + + Basis, BasisNotOrthonormal, TensorProductBasis, + basis_for_space, orthonormal_basis_for_space, monomial_basis_for_space) from modepy.nodes import ( equidistant_nodes, warp_and_blend_nodes, - tensor_product_nodes, legendre_gauss_lobatto_tensor_product_nodes) + tensor_product_nodes, legendre_gauss_lobatto_tensor_product_nodes, + + node_tuples_for_space, + equispaced_nodes_for_space, edge_clustered_nodes_for_space, + random_nodes_for_shape) from modepy.matrices import (vandermonde, resampling_matrix, differentiation_matrices, diff_matrix_permutation, inverse_mass_matrix, mass_matrix, - modal_face_mass_matrix, nodal_face_mass_matrix) -from modepy.quadrature import Quadrature, QuadratureRuleUnavailable + modal_face_mass_matrix, nodal_face_mass_matrix, + modal_mass_matrix_for_face, nodal_mass_matrix_for_face) +from modepy.quadrature import ( + Quadrature, QuadratureRuleUnavailable, + TensorProductQuadrature, LegendreGaussTensorProductQuadrature, + quadrature_for_space) from modepy.quadrature.jacobi_gauss import ( JacobiGaussQuadrature, LegendreGaussQuadrature, ChebyshevGaussQuadrature, - GaussGegenbauerQuadrature) + GaussGegenbauerQuadrature, + ) from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature @@ -56,6 +77,12 @@ __all__ = [ "__version__", + "Shape", "Face", "Simplex", "Hypercube", + "biunit_vertices_for_shape", "faces_for_shape", + "submesh_for_shape", + + "FunctionSpace", "PN", "QN", "space_for_shape", + "jacobi", "grad_jacobi", "simplex_onb", "grad_simplex_onb", "simplex_onb_with_mode_ids", "simplex_monomial_basis", "grad_simplex_monomial_basis", @@ -63,22 +90,34 @@ "simplex_best_available_basis", "grad_simplex_best_available_basis", "tensor_product_basis", "grad_tensor_product_basis", "legendre_tensor_product_basis", "grad_legendre_tensor_product_basis", + "symbolicize_function", + + "Basis", "BasisNotOrthonormal", "TensorProductBasis", + "basis_for_space", "orthonormal_basis_for_space", "monomial_basis_for_space", "equidistant_nodes", "warp_and_blend_nodes", "tensor_product_nodes", "legendre_gauss_lobatto_tensor_product_nodes", + "node_tuples_for_space", + "edge_clustered_nodes_for_space", "equispaced_nodes_for_space", + "random_nodes_for_shape", "vandermonde", "resampling_matrix", "differentiation_matrices", "diff_matrix_permutation", - "inverse_mass_matrix", "mass_matrix", "modal_face_mass_matrix", - "nodal_face_mass_matrix", + "inverse_mass_matrix", "mass_matrix", + "modal_face_mass_matrix", "nodal_face_mass_matrix", + "modal_mass_matrix_for_face", "nodal_mass_matrix_for_face", "Quadrature", "QuadratureRuleUnavailable", + "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature", + "quadrature_for_space", + "JacobiGaussQuadrature", "LegendreGaussQuadrature", "GaussLegendreQuadrature", "ChebyshevGaussQuadrature", "GaussGegenbauerQuadrature", "XiaoGimbutasSimplexQuadrature", "GrundmannMoellerSimplexQuadrature", "VioreanuRokhlinSimplexQuadrature", "ClenshawCurtisQuadrature", "FejerQuadrature", + "WitherdenVincentQuadrature", ] diff --git a/modepy/matrices.py b/modepy/matrices.py index 52d2fa1f..b237e7ce 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -21,9 +21,14 @@ """ +from warnings import warn import numpy as np import numpy.linalg as la +from modepy.shapes import Face, Simplex +from modepy.spaces import PN +from modepy.quadrature import Quadrature + __doc__ = r""" .. currentmodule:: modepy @@ -48,12 +53,9 @@ where :math:`(\phi_i)_i` is the basis of functions underlying :math:`V`. .. autofunction:: inverse_mass_matrix - .. autofunction:: mass_matrix - -.. autofunction:: modal_face_mass_matrix - -.. autofunction:: nodal_face_mass_matrix +.. autofunction:: modal_mass_matrix_for_face +.. autofunction:: nodal_mass_matrix_for_face Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives @@ -75,7 +77,7 @@ def vandermonde(functions, nodes): *functions* are allowed to return :class:`tuple` instances. In this case, a tuple of matrices is returned--i.e. this function - works directly on :func:`modepy.grad_simplex_onb` and returns + works directly on :func:`modepy.Basis.gradients` and returns a tuple of matrices. """ @@ -110,7 +112,7 @@ def resampling_matrix(basis, new_nodes, old_nodes, least_squares_ok=False): :arg basis: A sequence of basis functions accepting arrays of shape *(dims, npts)*, like those returned by - :func:`modepy.simplex_onb`. + :func:`modepy.orthonormal_basis_for_space`. :arg new_nodes: An array of shape *(dims, n_new_nodes)* :arg old_nodes: An array of shape *(dims, n_old_nodes)* :arg least_squares_ok: If *False*, then nodal values at *old_nodes* @@ -148,8 +150,9 @@ def is_square(m): np.dot(resample_vdm_new, la.pinv(vdm_old)), order="C") else: - raise RuntimeError("number of input nodes and number " - "of basis functions " + raise RuntimeError( + f"number of input nodes ({old_nodes.shape[1]}) " + f"and number of basis functions ({len(basis)}) " "do not agree--perhaps use least_squares_ok") @@ -160,10 +163,10 @@ def differentiation_matrices(basis, grad_basis, nodes, from_nodes=None): :arg basis: A sequence of basis functions accepting arrays of shape *(dims, npts)*, - like those returned by :func:`modepy.simplex_onb`. + like those returned by :func:`modepy.Basis.functions`. :arg grad_basis: A sequence of functions returning the gradients of *basis*, - like those returned by :func:`modepy.grad_simplex_onb`. + like those in :attr:`modepy.Basis.gradients`. :arg nodes: An array of shape *(dims, n_nodes)* :arg from_nodes: An array of shape *(dims, n_from_nodes)*. If *None*, assumed to be the same as *nodes*. @@ -242,44 +245,62 @@ def mass_matrix(basis, nodes): return la.inv(inverse_mass_matrix(basis, nodes)) -class _FaceMap: - def __init__(self, face_vertices): - """ - :arg face_vertices: an array of shape ``[dim, npts]``, where *npts* - should equal *dim*. - """ - vol_dim, npts = face_vertices.shape - if npts != vol_dim: - raise ValueError("face_vertices has wrong shape") +def modal_mass_matrix_for_face(face: Face, face_quad: Quadrature, + trial_functions, test_functions): + """ + .. versionadded:: 2020.3 + """ - self.origin = face_vertices[:, 0] - self.span = face_vertices[:, 1:] - self.origin[:, np.newaxis] + mapped_nodes = face.map_to_volume(face_quad.nodes) - self.face_dim = vol_dim - 1 + result = np.empty((len(test_functions), len(trial_functions))) - def __call__(self, points): - return (self.origin[:, np.newaxis] - + np.einsum("ad,dn->an", self.span, points*0.5 + 0.5)) + for i, test_f in enumerate(test_functions): + test_vals = test_f(mapped_nodes) + for j, trial_f in enumerate(trial_functions): + result[i, j] = (test_vals*trial_f(face_quad.nodes)) @ face_quad.weights + + return result +def nodal_mass_matrix_for_face(face: Face, face_quad: Quadrature, + trial_functions, test_functions, volume_nodes, face_nodes): + """ + .. versionadded :: 2020.3 + """ + face_vdm = vandermonde(trial_functions, face_nodes) + vol_vdm = vandermonde(test_functions, volume_nodes) + + modal_fmm = modal_mass_matrix_for_face( + face, face_quad, trial_functions, test_functions) + return la.inv(vol_vdm.T).dot(modal_fmm).dot(la.pinv(face_vdm)) + + +# {{{ deprecated junk + def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None): """ - :arg face_vertices: an array of shape ``[dim, npts]``, where *npts* - should equal *dim*. + :arg face_vertices: an array of shape ``(dims, nvertices)``. .. versionadded :: 2016.1 """ + warn("modal_face_mass_matrix is deprecated and will go away in 2022. " + "Use modal_mass_matrix_for_face instead.", + DeprecationWarning, stacklevel=2) + if test_basis is None: test_basis = trial_basis - fmap = _FaceMap(face_vertices) + vol_dims = face_vertices.shape[0] + + from modepy.quadrature import quadrature_for_space + quad = quadrature_for_space(PN(vol_dims - 1, order*2), Simplex(vol_dims-1)) - from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature - quad = GrundmannMoellerSimplexQuadrature(order, fmap.face_dim) - assert quad.exact_to > order*2 + assert quad.exact_to >= order*2 - mapped_nodes = fmap(quad.nodes) + from modepy.shapes import _simplex_face_to_vol_map + mapped_nodes = _simplex_face_to_vol_map(face_vertices, quad.nodes) nrows = len(test_basis) ncols = len(trial_basis) @@ -288,9 +309,7 @@ def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None): for i, test_f in enumerate(test_basis): test_vals = test_f(mapped_nodes) for j, trial_f in enumerate(trial_basis): - trial_vals = trial_f(mapped_nodes) - - result[i, j] = (test_vals*trial_vals).dot(quad.weights) + result[i, j] = (test_vals*trial_f(mapped_nodes)).dot(quad.weights) return result @@ -298,23 +317,28 @@ def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None): def nodal_face_mass_matrix(trial_basis, volume_nodes, face_nodes, order, face_vertices, test_basis=None): """ - :arg face_vertices: an array of shape ``[dim, npts]``, where *npts* - should equal *dim*. + :arg face_vertices: an array of shape ``(dims, nvertices)``. .. versionadded :: 2016.1 """ + warn("nodal_face_mass_matrix is deprecated and will go away in 2022. " + "Use nodal_mass_matrix_for_face instead.", + DeprecationWarning, stacklevel=2) + if test_basis is None: test_basis = trial_basis - fmap = _FaceMap(face_vertices) - - face_vdm = vandermonde(trial_basis, fmap(face_nodes)) # /!\ non-square + from modepy.shapes import _simplex_face_to_vol_map + face_vdm = vandermonde( + trial_basis, + _simplex_face_to_vol_map(face_vertices, face_nodes)) vol_vdm = vandermonde(test_basis, volume_nodes) modal_fmm = modal_face_mass_matrix( trial_basis, order, face_vertices, test_basis=test_basis) return la.inv(vol_vdm.T).dot(modal_fmm).dot(la.pinv(face_vdm)) +# }}} # vim: foldmethod=marker diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py index b255e6b4..8d456378 100644 --- a/modepy/modal_decay.py +++ b/modepy/modal_decay.py @@ -25,7 +25,7 @@ import numpy.linalg as la __doc__ = """Estimate the smoothness of a function represented in a basis -returned by :func:`modepy.simplex_onb`. +returned by :func:`modepy.orthonormal_basis_for_space`. The method implemented in this module follows this article: diff --git a/modepy/modes.py b/modepy/modes.py index f72d557d..52680048 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -22,13 +22,30 @@ """ +from warnings import warn import sys from math import sqrt +from functools import singledispatch, partial + import numpy as np +from modepy.spaces import FunctionSpace, PN, QN + + +__doc__ = """This functionality provides sets of basis functions for the +reference elements in :mod:`modepy.shapes`. + +.. currentmodule:: modepy + +Basis Retrieval +--------------- -__doc__ = """:mod:`modepy.modes` provides orthonormal bases and their -derivatives on unit simplices. +.. autoexception:: BasisNotOrthonormal +.. autoclass:: Basis + +.. autofunction:: basis_for_space +.. autofunction:: orthonormal_basis_for_space +.. autofunction:: monomial_basis_for_space Jacobi polynomials ------------------ @@ -36,11 +53,21 @@ .. currentmodule:: modepy .. autofunction:: jacobi(alpha, beta, n, x) - .. autofunction:: grad_jacobi(alpha, beta, n, x) -Dimension-independent basis getters for simplices -------------------------------------------------- +Conversion to Symbolic +---------------------- +.. autofunction:: symbolicize_function + +Tensor product adapter +---------------------- + +.. autoclass:: TensorProductBasis + +PKDO basis functions +-------------------- + +.. currentmodule:: modepy.modes .. |proriol-ref| replace:: Proriol, Joseph. "Sur une famille de polynomes á deux variables orthogonaux @@ -55,23 +82,6 @@ Scientific Computing 6, no. 4 (December 1, 1991): 345–390. http://dx.doi.org/10.1007/BF01060030 -.. autofunction:: simplex_onb_with_mode_ids -.. autofunction:: simplex_onb -.. autofunction:: grad_simplex_onb -.. autofunction:: simplex_monomial_basis_with_mode_ids -.. autofunction:: simplex_monomial_basis -.. autofunction:: grad_simplex_monomial_basis - -Dimension-independent basis getters for tensor-product bases ------------------------------------------------------------- - -.. autofunction:: tensor_product_basis -.. autofunction:: grad_tensor_product_basis - -Dimension-specific functions ----------------------------- - -.. currentmodule:: modepy.modes .. autofunction:: pkdo_2d .. autofunction:: grad_pkdo_2d @@ -84,9 +94,14 @@ .. autofunction:: monomial .. autofunction:: grad_monomial -Conversion to Symbolic ----------------------- -.. autofunction:: symbolicize_function +Redirections to Canonical Names +------------------------------- + +.. currentmodule:: modepy.modes + +.. class:: Basis + + See :class:`modepy.Basis`. """ @@ -251,10 +266,10 @@ def grad_pkdo_2d(order, rs): a, b = _rstoab(*rs) i, j = order - fa = jacobi(0, 0, i, a) - dfa = grad_jacobi(0, 0, i, a) - gb = jacobi(2*i+1, 0, j, b) - dgb = grad_jacobi(2*i+1, 0, j, b) + fa = _cse(jacobi(0, 0, i, a), f"leg_{i}") + dfa = _cse(grad_jacobi(0, 0, i, a), "dleg_{i}") + gb = _cse(jacobi(2*i+1, 0, j, b), f"jac_{2*i+1}_{j}") + dgb = _cse(grad_jacobi(2*i+1, 0, j, b), f"djac_{2*i+1}_{j}") # r-derivative # d/dr @@ -341,12 +356,12 @@ def grad_pkdo_3d(order, rst): a, b, c = _rsttoabc(*rst) i, j, k = order - fa = jacobi(0, 0, i, a) - dfa = grad_jacobi(0, 0, i, a) - gb = jacobi(2*i+1, 0, j, b) - dgb = grad_jacobi(2*i+1, 0, j, b) - hc = jacobi(2*(i+j)+2, 0, k, c) - dhc = grad_jacobi(2*(i+j)+2, 0, k, c) + fa = _cse(jacobi(0, 0, i, a), f"leg_{i}") + dfa = _cse(grad_jacobi(0, 0, i, a), f"dleg_{i}") + gb = _cse(jacobi(2*i+1, 0, j, b), f"jac_{2*i+1}") + dgb = _cse(grad_jacobi(2*i+1, 0, j, b), f"djac_{2*i+1}") + hc = _cse(jacobi(2*(i+j)+2, 0, k, c), f"jac_{2*(i+j)+2}") + dhc = _cse(grad_jacobi(2*(i+j)+2, 0, k, c), f"djac_{2*(i+j)+2}") # r-derivative # d/dr = da/dr d/da + db/dr d/db + dc/dr d/dx @@ -440,7 +455,7 @@ def diff_monomial(r, o): # }}} -# {{{ dimension-independent interface for simplices +# {{{ DEPRECATED dimension-independent interface for simplices def simplex_onb_with_mode_ids(dims, n): """Return a list of orthonormal basis functions in dimension *dims* of maximal @@ -460,39 +475,26 @@ def simplex_onb_with_mode_ids(dims, n): * |koornwinder-ref| * |dubiner-ref| - ... versionadded: 2018.1 + .. versionadded:: 2018.1 """ - from functools import partial - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - - if dims == 0: - def zerod_basis(x): - if len(x.shape) == 1: - return 1 - else: - return np.ones(x.shape[1]) + warn("simplex_onb_with_mode_ids is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) - return ((0,),), (zerod_basis,) - - elif dims == 1: + if dims == 1: mode_ids = tuple(range(n+1)) return mode_ids, tuple(partial(jacobi, 0, 0, i) for i in mode_ids) - elif dims == 2: - mode_ids = tuple(gnitstam(n, dims)) - return mode_ids, tuple(partial(pkdo_2d, order) for order in mode_ids) - elif dims == 3: - mode_ids = tuple(gnitstam(n, dims)) - return mode_ids, tuple(partial(pkdo_3d, order) for order in mode_ids) else: - raise NotImplementedError("%d-dimensional bases" % dims) + b = _SimplexONB(dims, n) + return b.mode_ids, b.functions def simplex_onb(dims, n): """Return a list of orthonormal basis functions in dimension *dims* of maximal total degree *n*. - :returns: a class:`tuple` of functions, each of which + :returns: a :class:`tuple` of functions, each of which accepts arrays of shape *(dims, npts)* and return the function values as an array of size *npts*. 'Scalar' evaluation, by passing just one vector of length *dims*, @@ -508,6 +510,11 @@ def simplex_onb(dims, n): Made return value a tuple, to make bases hashable. """ + warn("simplex_onb is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + mode_ids, basis = simplex_onb_with_mode_ids(dims, n) return basis @@ -532,7 +539,11 @@ def grad_simplex_onb(dims, n): Made return value a tuple, to make bases hashable. """ - from functools import partial + warn("grad_simplex_onb is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ as gnitstam @@ -560,10 +571,14 @@ def simplex_monomial_basis_with_mode_ids(dims, n): .. versionadded:: 2018.1 """ - from functools import partial - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - mode_ids = tuple(gnitstam(n, dims)) + warn("simplex_monomial_basis_with_mode_ids is deprecated. " + "Use monomial_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + from modepy.nodes import node_tuples_for_space + mode_ids = node_tuples_for_space(PN(dims, n)) + return mode_ids, tuple(partial(monomial, order) for order in mode_ids) @@ -571,7 +586,7 @@ def simplex_monomial_basis(dims, n): """Return a list of monomial basis functions in dimension *dims* of maximal total degree *n*. - :returns: a class:`tuple` of functions, each of which + :returns: a :class:`tuple` of functions, each of which accepts arrays of shape *(dims, npts)* and return the function values as an array of size *npts*. 'Scalar' evaluation, by passing just one vector of length *dims*, @@ -597,33 +612,41 @@ def grad_simplex_monomial_basis(dims, n): .. versionadded:: 2016.1 """ - from functools import partial + warn("grad_simplex_monomial_basis_with_mode_ids is deprecated. " + "Use monomial_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ as gnitstam return tuple(partial(grad_monomial, order) for order in gnitstam(n, dims)) -# undocumented for now def simplex_best_available_basis(dims, n): - if dims <= 3: - return simplex_onb(dims, n) - else: - return simplex_monomial_basis(dims, n) + warn("simplex_best_available_basis is deprecated. " + "Use basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + return basis_for_space(PN(dims, n)).functions -# undocumented for now def grad_simplex_best_available_basis(dims, n): - if dims <= 3: - return grad_simplex_onb(dims, n) - else: - return grad_simplex_monomial_basis(dims, n) + warn("grad_simplex_best_available_basis is deprecated. " + "Use basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + return basis_for_space(PN(dims, n)).gradients # }}} -# {{{ tensor product basis +# {{{ tensor product basis helpers class _TensorProductBasisFunction: + # multi_index is here just for debugging. + def __init__(self, multi_index, per_dim_functions): self.multi_index = multi_index self.per_dim_functions = per_dim_functions @@ -637,6 +660,8 @@ def __call__(self, x): class _TensorProductGradientBasisFunction: + # multi_index is here just for debugging. + def __init__(self, multi_index, per_dim_derivatives): self.multi_index = multi_index self.per_dim_derivatives = tuple(per_dim_derivatives) @@ -649,6 +674,10 @@ def __call__(self, x): return tuple(result) +# }}} + + +# {{{ DEPRECATED dimension-independent basis getters def tensor_product_basis(dims, basis_1d): """Adapt any iterable *basis_1d* of 1D basis functions into a *dims*-dimensional @@ -658,10 +687,17 @@ def tensor_product_basis(dims, basis_1d): .. versionadded:: 2017.1 """ - from pytools import generate_nonnegative_integer_tuples_below as gnitb + warn("tensor_product_basis is deprecated. " + "Use TensorProductBasis instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + from modepy.nodes import node_tuples_for_space + mode_ids = node_tuples_for_space(QN(dims, len(basis_1d))) + return tuple( _TensorProductBasisFunction(order, [basis_1d[i] for i in order]) - for order in gnitb(len(basis_1d), dims)) + for order in mode_ids) def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d): @@ -673,9 +709,14 @@ def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d): .. versionadded:: 2020.2 """ - from pytools import ( - wandering_element, - generate_nonnegative_integer_tuples_below as gnitb) + warn("grad_tensor_product_basis is deprecated. " + "Use TensorProductBasis instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + from pytools import wandering_element + from modepy.nodes import node_tuples_for_space + mode_ids = node_tuples_for_space(QN(dims, len(basis_1d))) func = (basis_1d, grad_basis_1d) return tuple( @@ -683,17 +724,25 @@ def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d): [func[i][k] for i, k in zip(iderivative, order)] for iderivative in wandering_element(dims) ]) - for order in gnitb(len(basis_1d), dims)) + for order in mode_ids) def legendre_tensor_product_basis(dims, order): - from functools import partial + warn("legendre_tensor_product_basis is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + basis = [partial(jacobi, 0, 0, n) for n in range(order + 1)] return tensor_product_basis(dims, basis) def grad_legendre_tensor_product_basis(dims, order): - from functools import partial + warn("grad_legendre_tensor_product_basis is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + basis = [partial(jacobi, 0, 0, n) for n in range(order + 1)] grad_basis = [partial(grad_jacobi, 0, 0, n) for n in range(order + 1)] return grad_tensor_product_basis(dims, basis, grad_basis) @@ -703,22 +752,22 @@ def grad_legendre_tensor_product_basis(dims, order): # {{{ conversion to symbolic -def symbolicize_function(f, dims, ref_coord_var_name="r"): +def symbolicize_function(f, dim, ref_coord_var_name="r"): """For a function *f* (basis or gradient) returned by one of the functions in this module, return a :mod:`pymbolic` expression representing the same function. - :arg dims: the number of dimensions of the reference element on which + :arg dim: the number of dimensions of the reference element on which *basis* is defined. .. versionadded:: 2020.2 """ import pymbolic.primitives as p - r_sym = p.make_sym_vector(ref_coord_var_name, dims) + r_sym = p.make_sym_vector(ref_coord_var_name, dim) result = f(r_sym) - if dims == 1: + if dim == 1: # Work around inconsistent 1D stupidity. Grrrr! # (We fed it an object array, and it gave one back, i.e. it treated its # argument as a scalar instead of indexing into it. That tends to @@ -735,4 +784,271 @@ def symbolicize_function(f, dims, ref_coord_var_name="r"): # }}} + +# {{{ basis interface + +class BasisNotOrthonormal(Exception): + pass + + +class Basis: + """ + .. automethod:: orthonormality_weight(r) + .. autoattribute:: mode_ids + .. autoattribute:: functions + .. autoattribute:: gradients + """ + + def orthonormality_weight(self, rvec): + raise NotImplementedError + + @property + def mode_ids(self): + """Return a tuple of of mode (basis function) identifiers, one for + each basis function. Mode identifiers should generally be viewed as opaque. + They are hashable. For some bases, they are tuples of length matching + the number of dimensions and indicating the order along each reference + axis. + """ + raise NotImplementedError + + @property + def functions(self): + """Return a tuple of (callable) basis functions of length matching + :attr:`mode_ids`. Each function takes a vector of (r,s,t) reference + coordinates as input. + """ + raise NotImplementedError + + @property + def gradients(self): + """Return a tuple of basis functions of length matching :attr:`mode_ids`. + Each function takes a vector of (r,s,t) reference coordinates as input. + """ + raise NotImplementedError + +# }}} + + +# {{{ space-based basis retrieval + +@singledispatch +def basis_for_space(space: FunctionSpace) -> Basis: + raise NotImplementedError(type(space).__name__) + + +@singledispatch +def orthonormal_basis_for_space(space: FunctionSpace) -> Basis: + raise NotImplementedError(type(space).__name__) + + +@singledispatch +def monomial_basis_for_space(space: FunctionSpace) -> Basis: + raise NotImplementedError(type(space).__name__) + +# }}} + + +def zerod_basis(x): + assert len(x) == 0 + x_sub = np.ones(x.shape[1:], x.dtype) + return 1 + x_sub + + +# {{{ PN bases + +def _pkdo_1d(order, r): + i, = order + r0, = r + return jacobi(0, 0, i, r0) + + +def _grad_pkdo_1d(order, r): + i, = order + r0, = r + return (grad_jacobi(0, 0, i, r0),) + + +class _SimplexBasis(Basis): + def __init__(self, dim, order): + self._dim = dim + self._order = order + + assert isinstance(dim, int) + assert isinstance(order, int) + + @property + def mode_ids(self): + from pytools import \ + generate_nonnegative_integer_tuples_summing_to_at_most as gnitsam + return tuple(gnitsam(self._order, self._dim)) + + +class _SimplexONB(_SimplexBasis): + is_orthonormal = True + + def orthonormality_weight(self, rvec): + return 1 + + @property + def functions(self): + if self._dim == 0: + return (zerod_basis,) + elif self._dim == 1: + return tuple(partial(_pkdo_1d, mid) for mid in self.mode_ids) + elif self._dim == 2: + return tuple(partial(pkdo_2d, mid) for mid in self.mode_ids) + elif self._dim == 3: + return tuple(partial(pkdo_3d, mid) for mid in self.mode_ids) + else: + raise NotImplementedError(f"basis in {self._dim} dimensions") + + @property + def gradients(self): + if self._dim == 1: + return tuple(partial(_grad_pkdo_1d, mid) for mid in self.mode_ids) + elif self._dim == 2: + return tuple(partial(grad_pkdo_2d, mid) for mid in self.mode_ids) + elif self._dim == 3: + return tuple(partial(grad_pkdo_3d, mid) for mid in self.mode_ids) + else: + raise NotImplementedError(f"gradient in {self._dim} dimensions") + + +class _SimplexMonomialBasis(_SimplexBasis): + def orthonormality_weight(self, rvec): + raise BasisNotOrthonormal + + @property + def functions(self): + return tuple(partial(monomial, mid) for mid in self.mode_ids) + + @property + def gradients(self): + return tuple(partial(grad_monomial, mid) for mid in self.mode_ids) + + +@basis_for_space.register(PN) +def _(space: PN): + if space.spatial_dim <= 3: + return _SimplexONB(space.spatial_dim, space.order) + else: + return _SimplexMonomialBasis(space.spatial_dim, space.order) + + +@orthonormal_basis_for_space.register(PN) +def _(space: PN): + return _SimplexONB(space.spatial_dim, space.order) + + +@monomial_basis_for_space.register(PN) +def _(space: PN): + return _SimplexMonomialBasis(space.spatial_dim, space.order) + +# }}} + + +# {{{ QN bases + +class TensorProductBasis(Basis): + """Adapts multiple one-dimensional bases into a tensor product basis. + + .. automethod:: __init__ + """ + + def __init__(self, bases_1d, grad_bases_1d, orth_weight): + """ + :arg bases_1d: a sequence (one entry per axis/dimension) + of sequences (representing the basis) of 1D functions + representing the approximation basis. + :arg grad_bases_1d: a sequence (one entry per axis/dimension) + representing the derivatives of *bases_1d*. + """ + self._bases_1d = bases_1d + self._grad_bases_1d = grad_bases_1d + self._orth_weight = orth_weight + + if len(bases_1d) != len(grad_bases_1d): + raise ValueError("bases_1d and grad_bases_1d must have the same length") + + for i, (b, gb) in enumerate(zip(bases_1d, grad_bases_1d)): + if len(b) != len(gb): + raise ValueError( + f"bases_1d[{i}] and grad_bases_1d[{i}] " + "must have the same length") + + def orthonormality_weight(self): + if self._orth_weight is None: + raise BasisNotOrthonormal + else: + return self._orth_weight + + @property + def _dim(self): + return len(self._bases_1d) + + @property + def mode_ids(self): + from pytools import generate_nonnegative_integer_tuples_below as gnitb + return tuple(gnitb([len(b) for b in self._bases_1d])) + + @property + def functions(self): + return tuple( + _TensorProductBasisFunction(mid, + [self._bases_1d[iaxis][mid_i] + for iaxis, mid_i in enumerate(mid)]) + for mid in self.mode_ids) + + @property + def gradients(self): + from pytools import wandering_element + func = (self._bases_1d, self._grad_bases_1d) + return tuple( + _TensorProductGradientBasisFunction(mid, [ + [func[is_deriv][iaxis][mid_i] + for iaxis, (is_deriv, mid_i) in + enumerate(zip(iderivative, mid))] + for iderivative in wandering_element(self._dim) + ]) + for mid in self.mode_ids) + + +@orthonormal_basis_for_space.register(QN) +def _(space: QN): + order = space.order + dim = space.spatial_dim + return TensorProductBasis( + [[partial(jacobi, 0, 0, n) for n in range(order + 1)]] * dim, + [[partial(grad_jacobi, 0, 0, n) for n in range(order + 1)]] * dim, + orth_weight=1) + + +@basis_for_space.register(QN) +def _(space: QN): + return orthonormal_basis_for_space(space) + + +def _monomial_1d(order, r): + return r**order + + +def _grad_monomial_1d(order, r): + if order == 0: + return 0*r + else: + return order*r**(order-1) + + +@monomial_basis_for_space.register(QN) +def _(space: QN): + order = space.order + dim = space.spatial_dim + return TensorProductBasis( + [[partial(_monomial_1d, n) for n in range(order + 1)]] * dim, + [[partial(_grad_monomial_1d, n) for n in range(order + 1)]] * dim, + orth_weight=None) + +# }}} + # vim: foldmethod=marker diff --git a/modepy/nodes.py b/modepy/nodes.py index e8d0dfcc..a2fc8146 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -1,3 +1,34 @@ +# {{{ docstring + +""" +Generic Shape-Based Interface +----------------------------- + +.. currentmodule:: modepy + +.. autofunction:: node_tuples_for_space +.. autofunction:: equispaced_nodes_for_space +.. autofunction:: edge_clustered_nodes_for_space +.. autofunction:: random_nodes_for_shape + +Simplices +--------- + +.. autofunction:: equidistant_nodes +.. autofunction:: warp_and_blend_nodes + +Also see :class:`modepy.VioreanuRokhlinSimplexQuadrature` if nodes on the +boundary are not required. + +Hypercubes +---------- + +.. currentmodule:: modepy + +.. autofunction:: tensor_product_nodes +.. autofunction:: legendre_gauss_lobatto_tensor_product_nodes +""" + __copyright__ = "Copyright (C) 2009, 2010, 2013 Andreas Kloeckner, " \ "Tim Warburton, Jan Hesthaven, Xueyu Zhu" @@ -21,9 +52,19 @@ THE SOFTWARE. """ +# }}} + +from typing import List, Tuple import numpy as np import numpy.linalg as la +from functools import singledispatch, partial + +from modepy.shapes import Shape, Simplex, Hypercube +from modepy.spaces import FunctionSpace, PN, QN + + +# {{{ equidistant nodes def equidistant_nodes(dims, n, node_tuples=None): """ @@ -33,30 +74,24 @@ def equidistant_nodes(dims, n, node_tuples=None): :arg node_tuples: a list of tuples of integers indicating the node order. Use default order if *None*, see :func:`pytools.generate_nonnegative_integer_tuples_summing_to_at_most`. - :returns: An array of shape *(dims, nnodes)* containing unit coordinates + :returns: An array of shape *(dims, nnodes)* containing bi-unit coordinates of the interpolation nodes. (see :ref:`tri-coords` and :ref:`tet-coords`) """ + space = PN(dims, n) if node_tuples is None: - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - node_tuples = list(gnitstam(n, dims)) + node_tuples = node_tuples_for_space(space) else: - try: - from math import comb # comb is v3.8+ - node_count = comb(n + dims, dims) - except ImportError: - from functools import reduce - from operator import mul - node_count = reduce(mul, range(n + 1, n + dims + 1), 1) \ - // reduce(mul, range(1, dims + 1), 1) - - if len(node_tuples) != node_count: + if len(node_tuples) != space.space_dim: raise ValueError("'node_tuples' list does not have the correct length") - # shape: (2, nnodes) + # shape: (dims, nnodes) return (np.array(node_tuples, dtype=np.float64)/n*2 - 1).T +# }}} + + +# {{{ warp and blend simplex nodes def warp_factor(n, output_nodes, scaled=True): """Compute warp function at order *n* and evaluate it at @@ -69,9 +104,9 @@ def warp_factor(n, output_nodes, scaled=True): equi_nodes = np.linspace(-1, 1, n+1) from modepy.matrices import vandermonde - from modepy.modes import simplex_onb + from modepy.modes import jacobi - basis = simplex_onb(1, n) + basis = [partial(jacobi, 0, 0, n) for n in range(n + 1)] Veq = vandermonde(basis, equi_nodes) # noqa # create interpolator from equi_nodes to output_nodes @@ -126,13 +161,12 @@ def warp_and_blend_nodes_2d(n, node_tuples=None): except IndexError: alpha = 5/3 + space = PN(2, n) if node_tuples is None: - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - node_tuples = list(gnitstam(n, 2)) + node_tuples = node_tuples_for_space(space) else: - if len(node_tuples) != (n+1)*(n+2)//2: - raise ValueError("node_tuples list does not have the correct length") + if len(node_tuples) != space.space_dim: + raise ValueError("'node_tuples' list does not have the correct length") # shape: (2, nnodes) unit_nodes = (np.array(node_tuples, dtype=np.float64)/n*2 - 1).T @@ -163,13 +197,12 @@ def warp_and_blend_nodes_3d(n, node_tuples=None): except IndexError: alpha = 1. + space = PN(3, n) if node_tuples is None: - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - node_tuples = list(gnitstam(n, 3)) + node_tuples = node_tuples_for_space(space) else: - if len(node_tuples) != (n+1)*(n+2)*(n+3)//6: - raise ValueError("node_tuples list does not have the correct length") + if len(node_tuples) != space.space_dim: + raise ValueError("'node_tuples' list does not have the correct length") # shape: (3, nnodes) unit_nodes = (np.array(node_tuples, dtype=np.float64)/n*2 - 1).T @@ -291,31 +324,152 @@ def warp_and_blend_nodes(dims, n, node_tuples=None): # }}} +# }}} + # {{{ tensor product nodes def tensor_product_nodes(dims, nodes_1d): """ - :returns: an array of shape ``(dims, nnodes_1d**dims)``. The - order of nodes is such that the nodes along the last - axis vary fastest. + :returns: an array of shape ``(dims, nnodes_1d**dims)``. .. versionadded:: 2017.1 + + .. versionchanged:: 2020.3 + + The node ordering has changed and is no longer documented. """ nnodes_1d = len(nodes_1d) result = np.empty((dims,) + (nnodes_1d,) * dims) for d in range(dims): - result[-d - 1] = nodes_1d.reshape(*((-1,) + (1,)*d)) + result[d] = nodes_1d.reshape(*((-1,) + (1,)*d)) return result.reshape(dims, -1) def legendre_gauss_lobatto_tensor_product_nodes(dims, n): from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes - return tensor_product_nodes(dims, - legendre_gauss_lobatto_nodes(n)) + return tensor_product_nodes(dims, legendre_gauss_lobatto_nodes(n)) + +# }}} + + +# {{{ space-based interface + +@singledispatch +def node_tuples_for_space(space: FunctionSpace) -> List[Tuple[int]]: + raise NotImplementedError(type(space).__name__) + + +@singledispatch +def equispaced_nodes_for_space(space: FunctionSpace, shape: Shape): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + + +@singledispatch +def edge_clustered_nodes_for_space(space: FunctionSpace, shape: Shape): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + + +@singledispatch +def random_nodes_for_shape(shape: Shape, nnodes: int, rng=None): + """ + :arg generator: a :class:`numpy.random.Generator`. + :returns: a :class:`numpy.ndarray` that returns an array of + shape `(dim, nnodes)` of random nodes in the reference element. + """ + raise NotImplementedError(type(shape).__name__) + +# }}} + + +# {{{ PN + +@node_tuples_for_space.register(PN) +def _(space: PN): + from pytools import \ + generate_nonnegative_integer_tuples_summing_to_at_most as gnitsam + return tuple(gnitsam(space.order, space.spatial_dim)) + + +@equispaced_nodes_for_space.register(PN) +def _(space: PN, shape: Simplex): + if not isinstance(shape, Simplex): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + return (np.array(node_tuples_for_space(space), dtype=np.float64) + / space.order*2 - 1).T + + +@edge_clustered_nodes_for_space.register(PN) +def _(space: PN, shape: Simplex): + if not isinstance(shape, Simplex): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + return warp_and_blend_nodes(space.spatial_dim, space.order) + + +@random_nodes_for_shape.register(Simplex) +def _(shape: Simplex, nnodes: int, rng=None): + if rng is None: + rng = np.random + + result = np.zeros((shape.dim, nnodes)) + nnodes_obtained = 0 + while True: + new_nodes = rng.uniform(-1.0, 1.0, size=(shape.dim, nnodes-nnodes_obtained)) + new_nodes = new_nodes[:, new_nodes.sum(axis=0) < 2-shape.dim] + nnew_nodes = new_nodes.shape[1] + result[:, nnodes_obtained:nnodes_obtained+nnew_nodes] = new_nodes + nnodes_obtained += nnew_nodes + + if nnodes_obtained == nnodes: + return result # }}} +# {{{ QN + +@node_tuples_for_space.register(QN) +def _(space: QN): + from pytools import \ + generate_nonnegative_integer_tuples_below as gnitb + return tuple(gnitb(space.order, space.spatial_dim)) + + +@equispaced_nodes_for_space.register(QN) +def _(space: QN, shape: Hypercube): + if not isinstance(shape, Hypercube): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + return (np.array(node_tuples_for_space(space), dtype=np.float64) + / space.order*2 - 1).T + + +@edge_clustered_nodes_for_space.register(QN) +def _(space: QN, shape: Hypercube): + if not isinstance(shape, Hypercube): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + return legendre_gauss_lobatto_tensor_product_nodes( + space.spatial_dim, space.order) + + +@random_nodes_for_shape.register(Hypercube) +def _(shape: Hypercube, nnodes: int, rng=None): + if rng is None: + rng = np.random + return rng.uniform(-1.0, 1.0, size=(shape.dim, nnodes)) + +# }}} + # vim: foldmethod=marker diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 7ca637c8..67f37be9 100644 --- a/modepy/quadrature/__init__.py +++ b/modepy/quadrature/__init__.py @@ -1,3 +1,20 @@ +""" +.. currentmodule:: modepy + +.. autoclass:: Quadrature + +.. autofunction:: quadrature_for_space + +Redirections to Canonical Names +------------------------------- + +.. currentmodule:: modepy.quadrature + +.. class:: Quadrature + + See :class:`modepy.Quadrature`. +""" + __copyright__ = ("Copyright (C) 2009, 2010, 2013 Andreas Kloeckner, Tim Warburton, " "Jan Hesthaven, Xueyu Zhu") @@ -21,8 +38,11 @@ THE SOFTWARE. """ +from functools import singledispatch import numpy as np +from modepy.shapes import Shape, Simplex, Hypercube +from modepy.spaces import FunctionSpace, PN, QN class QuadratureRuleUnavailable(RuntimeError): @@ -82,3 +102,84 @@ def __init__(self, quad, left, right): Quadrature.__init__(self, left + (quad.nodes+1) / 2 * length, quad.weights * half_length) + + +class TensorProductQuadrature(Quadrature): + """ + .. automethod:: __init__ + """ + + def __init__(self, dims, quad): + """ + :arg quad: a :class:`Quadrature` class for one-dimensional intervals. + """ + + from modepy.nodes import tensor_product_nodes + x = tensor_product_nodes(dims, quad.nodes) + from itertools import product + w = np.fromiter( + (np.prod(w) for w in product(quad.weights, repeat=dims)), + dtype=np.float, + count=quad.weights.size**dims) + assert w.size == x.shape[1] + + super().__init__(x, w) + self.exact_to = quad.exact_to + + +class LegendreGaussTensorProductQuadrature(TensorProductQuadrature): + def __init__(self, N, dims, backend=None): # noqa: N803 + from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature + super().__init__( + dims, LegendreGaussQuadrature(N, backend=backend)) + + +# {{{ quadrature + +@singledispatch +def quadrature_for_space(space: FunctionSpace, shape: Shape) -> Quadrature: + """ + :returns: a :class:`~modepy.Quadrature` that exactly integrates the functions + in *space*. + """ + raise NotImplementedError((type(space).__name__, type(shape).__name)) + + +@quadrature_for_space.register(PN) +def _(space: PN, shape: Simplex): + if not isinstance(shape, Simplex): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + import modepy as mp + try: + quad = mp.XiaoGimbutasSimplexQuadrature(space.order, space.spatial_dim) + except mp.QuadratureRuleUnavailable: + quad = mp.GrundmannMoellerSimplexQuadrature( + space.order//2, space.spatial_dim) + + assert quad.exact_to >= space.order + + return quad + + +@quadrature_for_space.register(QN) +def _(space: QN, shape: Hypercube): + if not isinstance(shape, Hypercube): + raise NotImplementedError((type(space).__name__, type(shape).__name)) + if space.spatial_dim != shape.dim: + raise ValueError("spatial dimensions of shape and space must match") + + import modepy as mp + if space.spatial_dim == 0: + quad = mp.Quadrature(np.empty((0, 1)), np.empty((0, 1))) + else: + from modepy.quadrature import LegendreGaussTensorProductQuadrature + quad = LegendreGaussTensorProductQuadrature(space.order, space.spatial_dim) + + return quad + +# }}} + +# vim: foldmethod=marker diff --git a/modepy/quadrature/witherden_vincent.py b/modepy/quadrature/witherden_vincent.py index 9d8b9a4c..2731d083 100644 --- a/modepy/quadrature/witherden_vincent.py +++ b/modepy/quadrature/witherden_vincent.py @@ -28,15 +28,15 @@ class WitherdenVincentQuadrature(Quadrature): hexahedra. The integration domain is the unit hypercube :math:`[-1, 1]^d`, where :math:`d` - is the dimension. The quadrature rules are adapted from:: + is the dimension. The quadrature rules are adapted from: F. D. Witherden, P. E. Vincent, On the Identification of Symmetric Quadrature Rules for Finite Element Methods, Computers & Mathematics with Applications, Vol. 69, pp. 1232--1241, 2015, - `10.1016/j.camwa.2015.03.017 http://dx.doi.org/10.1016/j.camwa.2015.03.017`_. + `DOI `_. - .. versionadded: 2020.5 + .. versionadded: 2020.3 .. automethod:: __init__ .. automethod:: __call__ diff --git a/modepy/quadrature/xiao_gimbutas.py b/modepy/quadrature/xiao_gimbutas.py index 31e7bbe7..caa98e58 100644 --- a/modepy/quadrature/xiao_gimbutas.py +++ b/modepy/quadrature/xiao_gimbutas.py @@ -64,7 +64,7 @@ def __init__(self, order, dims): elif dims == 3: from modepy.quadrature.xg_quad_data import tetrahedron_table as table else: - raise ValueError("invalid dimensionality") + raise QuadratureRuleUnavailable("invalid dimensionality") try: order_table = table[order] except KeyError: diff --git a/modepy/shapes.py b/modepy/shapes.py new file mode 100644 index 00000000..5cf162b8 --- /dev/null +++ b/modepy/shapes.py @@ -0,0 +1,547 @@ +# {{{ docstring + +r""" +:mod:`modepy.shapes` provides a generic description of the supported shapes +(i.e. reference elements). + +.. currentmodule:: modepy + +.. autoclass:: Shape +.. autoclass:: Face +.. autofunction:: biunit_vertices_for_shape +.. autofunction:: faces_for_shape + +Simplices +^^^^^^^^^ + +.. autoclass:: Simplex + +.. _tri-coords: + +Coordinates on the triangle +--------------------------- + +Bi-unit coordinates :math:`(r, s)` (also called 'unit' coordinates):: + + ^ s + | + C + |\ + | \ + | O + | \ + | \ + A-----B--> r + +Vertices in bi-unit coordinates:: + + O = ( 0, 0) + A = (-1, -1) + B = ( 1, -1) + C = (-1, 1) + +Equilateral coordinates :math:`(x, y)`:: + + C + / \ + / \ + / \ + / O \ + / \ + A-----------B + +Vertices in equilateral coordinates:: + + O = ( 0, 0) + A = (-1, -1/sqrt(3)) + B = ( 1, -1/sqrt(3)) + C = ( 0, 2/sqrt(3)) + +.. _tet-coords: + +Coordinates on the tetrahedron +------------------------------ + +Bi-unit coordinates :math:`(r, s, t)` (also called 'unit' coordinates):: + + ^ s + | + C + /|\ + / | \ + / | \ + / | \ + / O| \ + / __A-----B---> r + /_--^ ___--^^ + ,D--^^^ + t L + +(squint, and it might start making sense...) + +Vertices in bi-unit coordinates :math:`(r, s, t)`:: + + O = ( 0, 0, 0) + A = (-1, -1, -1) + B = ( 1, -1, -1) + C = (-1, 1, -1) + D = (-1, -1, 1) + +Vertices in equilateral coordinates :math:`(x, y, z)`:: + + O = ( 0, 0, 0) + A = (-1, -1/sqrt(3), -1/sqrt(6)) + B = ( 1, -1/sqrt(3), -1/sqrt(6)) + C = ( 0, 2/sqrt(3), -1/sqrt(6)) + D = ( 0, 0, 3/sqrt(6)) + +Hypercubes +^^^^^^^^^^ + +.. autoclass:: Hypercube + +.. _square-coords: + +Coordinates on the square +------------------------- + +Bi-unit coordinates on :math:`(r, s)` (also called 'unit' coordinates):: + + ^ s + | + C---------D + | | + | | + | O | + | | + | | + A---------B --> r + + +Vertices in bi-unit coordinates:: + + O = ( 0, 0) + A = (-1, -1) + B = ( 1, -1) + C = (-1, 1) + D = ( 1, 1) + +.. _cube-coords: + +Coordinates on the cube +----------------------- + +Bi-unit coordinates on :math:`(r, s, t)` (also called 'unit' coordinates):: + + t + ^ + | + E----------G + |\ |\ + | \ | \ + | \ | \ + | F------+---H + | | O | | + A---+------C---|--> s + \ | \ | + \ | \ | + \| \| + B----------D + \ + v r + +Vertices in bi-unit coordinates:: + + O = ( 0, 0, 0) + A = (-1, -1, -1) + B = ( 1, -1, -1) + C = (-1, 1, -1) + D = ( 1, 1, -1) + E = (-1, -1, 1) + F = ( 1, -1, 1) + G = (-1, 1, 1) + H = ( 1, 1, 1) + +The order of the vertices in the hypercubes follows binary counting +in ``tsr`` (i.e. in reverse axis order). +For example, in 3D, ``A, B, C, D, ...`` is ``000, 001, 010, 011, ...``. + +Submeshes +--------- +.. autofunction:: submesh_for_shape + +Redirections to Canonical Names +------------------------------- + +.. currentmodule:: modepy.shapes + +.. class:: Shape + + See :class:`modepy.Shape`. + +.. class:: Face + + See :class:`modepy.Face`. + +.. class:: Simplex + + See :class:`modepy.Simplex`. + +.. class:: Hypercube + + See :class:`modepy.Hypercube`. +""" + +# }}} + +__copyright__ = """ +Copyright (c) 2013 Andreas Kloeckner +Copyright (c) 2020 Alexandru Fikl +""" + +__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. +""" + +import numpy as np +from typing import Tuple, Callable + +from functools import singledispatch, partial +from dataclasses import dataclass + + +# {{{ interface + +@dataclass(frozen=True) +class Shape: + """ + .. attribute :: dim + .. attribute :: nfaces + .. attribute :: nvertices + """ + dim: int + + +@singledispatch +def biunit_vertices_for_shape(shape: Shape): + """ + :returns: an :class:`~numpy.ndarray` of shape `(dim, nvertices)`. + """ + raise NotImplementedError(type(shape).__name__) + + +@dataclass(frozen=True) +class Face: + """Mix-in to be used with a concrete :class:`Shape` subclass to represent + geometry information about a face of a shape. + + .. attribute:: volume_shape + + The volume :class:`Shape` from which this face descends. + + .. attribute:: face_index + + The face index in :attr:`volume_shape` of this face. + + .. attribute:: volume_vertex_indices + + A tuple of indices into the vertices returned by + :func:`biunit_vertices_for_shape` for the :attr:`volume_shape`. + + .. attribute:: map_to_volume + + A :class:`~collections.abc.Callable` that takes an array of + size `(dim, nnodes)` of unit nodes on the face represented by + *face_vertices* and maps them to the :attr:`volume_shape`. + """ + volume_shape: Shape + face_index: int + volume_vertex_indices: Tuple[int] + map_to_volume: Callable[[np.ndarray], np.ndarray] + + +@singledispatch +def faces_for_shape(shape: Shape): + """ + :results: a tuple of :class:`Face` representing the faces of *shape*. + """ + raise NotImplementedError(type(shape).__name__) + +# }}} + + +# {{{ simplex + +class Simplex(Shape): + @property + def nfaces(self): + return self.dim + 1 + + @property + def nvertices(self): + return self.dim + 1 + + +@dataclass(frozen=True) +class _SimplexFace(Simplex, Face): + pass + + +@biunit_vertices_for_shape.register(Simplex) +def _(shape: Simplex): + result = np.empty((shape.dim, shape.dim+1), np.float64) + result.fill(-1) + + for i in range(shape.dim): + result[i, i+1] = 1 + + return result + + +def _simplex_face_to_vol_map(face_vertices, p: np.ndarray): + dim, npoints = face_vertices.shape + if npoints != dim: + raise ValueError("'face_vertices' has wrong shape") + + origin = face_vertices[:, 0].reshape(-1, 1) + face_basis = face_vertices[:, 1:] - origin + + return origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2) + + +@faces_for_shape.register(Simplex) +def _(shape: Simplex): + face_vertex_indices = np.empty((shape.dim + 1, shape.dim), dtype=np.int) + indices = np.arange(shape.dim + 1) + + for iface in range(shape.nfaces): + face_vertex_indices[iface, :] = \ + np.hstack([indices[:iface], indices[iface + 1:]]) + + vertices = biunit_vertices_for_shape(shape) + return [ + _SimplexFace( + dim=shape.dim-1, + volume_shape=shape, face_index=iface, + volume_vertex_indices=tuple(fvi), + map_to_volume=partial(_simplex_face_to_vol_map, vertices[:, fvi])) + for iface, fvi in enumerate(face_vertex_indices)] + +# }}} + + +# {{{ hypercube + +class Hypercube(Shape): + @property + def nfaces(self): + return 2 * self.dim + + @property + def nvertices(self): + return 2**self.dim + + +@dataclass(frozen=True) +class _HypercubeFace(Hypercube, Face): + pass + + +@biunit_vertices_for_shape.register(Hypercube) +def _(shape: Hypercube): + from modepy.nodes import tensor_product_nodes + return tensor_product_nodes(shape.dim, np.array([-1.0, 1.0])) + + +def _hypercube_face_to_vol_map(face_vertices: np.ndarray, p: np.ndarray): + dim, npoints = face_vertices.shape + if npoints != 2**(dim - 1): + raise ValueError("'face_vertices' has wrong shape") + + origin = face_vertices[:, 0].reshape(-1, 1) + + # works up to (and including) 3D: + # - no-op for 1D, 2D + # - For square faces, eliminate middle node + face_basis = face_vertices[:, 1:3] - origin + + return origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2) + + +@faces_for_shape.register(Hypercube) +def _(shape: Hypercube): + # FIXME: replace by nicer n-dimensional formula + face_vertex_indices = { + 1: ((0b0,), (0b1,)), + 2: ((0b00, 0b01), (0b10, 0b11), (0b00, 0b10), (0b01, 0b11)), + 3: ( + (0b000, 0b001, 0b010, 0b011,), + (0b100, 0b101, 0b110, 0b111,), + + (0b000, 0b010, 0b100, 0b110,), + (0b001, 0b011, 0b101, 0b111,), + + (0b000, 0b001, 0b100, 0b101,), + (0b010, 0b011, 0b110, 0b111,), + ) + }[shape.dim] + + vertices = biunit_vertices_for_shape(shape) + return [ + _HypercubeFace( + dim=shape.dim-1, + volume_shape=shape, face_index=iface, + volume_vertex_indices=fvi, + map_to_volume=partial(_hypercube_face_to_vol_map, vertices[:, fvi])) + for iface, fvi in enumerate(face_vertex_indices)] + +# }}} + + +# {{{ submeshes + +@singledispatch +def submesh_for_shape(shape: Shape, node_tuples): + """Return a list of tuples of indices into the node list that + generate a tesselation of the reference element. + + :arg node_tuples: A list of tuples *(i, j, ...)* of integers + indicating node positions inside the unit element. The + returned list references indices in this list. + + :func:`modepy.node_tuples_for_space` may be used to generate *node_tuples*. + + .. versionadded:: 2020.3 + """ + raise NotImplementedError(type(shape).__name__) + + +@submesh_for_shape.register(Simplex) +def _(shape: Simplex, node_tuples): + from pytools import single_valued, add_tuples + dims = single_valued(len(nt) for nt in node_tuples) + + node_dict = { + ituple: idx + for idx, ituple in enumerate(node_tuples)} + + if dims == 1: + result = [] + + def try_add_line(d1, d2): + try: + result.append(( + node_dict[add_tuples(current, d1)], + node_dict[add_tuples(current, d2)], + )) + except KeyError: + pass + + for current in node_tuples: + try_add_line((0,), (1,),) + + return result + elif dims == 2: + # {{{ triangle sub-mesh + result = [] + + def try_add_tri(d1, d2, d3): + try: + result.append(( + node_dict[add_tuples(current, d1)], + node_dict[add_tuples(current, d2)], + node_dict[add_tuples(current, d3)], + )) + except KeyError: + pass + + for current in node_tuples: + # this is a tesselation of a square into two triangles. + # subtriangles that fall outside of the master triangle are + # simply not added. + + # positively oriented + try_add_tri((0, 0), (1, 0), (0, 1)) + try_add_tri((1, 0), (1, 1), (0, 1)) + + return result + + # }}} + elif dims == 3: + # {{{ tet sub-mesh + + def try_add_tet(d1, d2, d3, d4): + try: + result.append(( + node_dict[add_tuples(current, d1)], + node_dict[add_tuples(current, d2)], + node_dict[add_tuples(current, d3)], + node_dict[add_tuples(current, d4)], + )) + except KeyError: + pass + + result = [] + for current in node_tuples: + # this is a tesselation of a cube into six tets. + # subtets that fall outside of the master tet are simply not added. + + # positively oriented + try_add_tet((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) + try_add_tet((1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0)) + try_add_tet((1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1)) + + try_add_tet((1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0)) + try_add_tet((0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1)) + try_add_tet((0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0)) + + return result + + # }}} + else: + raise NotImplementedError("%d-dimensional sub-meshes" % dims) + + +@submesh_for_shape.register(Hypercube) +def _(shape: Hypercube, node_tuples): + from pytools import single_valued, add_tuples + dims = single_valued(len(nt) for nt in node_tuples) + + node_dict = { + ituple: idx + for idx, ituple in enumerate(node_tuples)} + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + result = [] + for current in node_tuples: + try: + result.append(tuple( + node_dict[add_tuples(current, offset)] + for offset in gnitb(2, dims))) + + except KeyError: + pass + + return result + +# }}} + + +# vim: foldmethod=marker diff --git a/modepy/spaces.py b/modepy/spaces.py new file mode 100644 index 00000000..f0898242 --- /dev/null +++ b/modepy/spaces.py @@ -0,0 +1,146 @@ +""" +.. currentmodule:: modepy + +Function Spaces +--------------- + +.. autoclass:: FunctionSpace +.. autoclass:: PN +.. autoclass:: QN +.. autoclass:: space_for_shape + +Redirections to Canonical Names +------------------------------- + +.. currentmodule:: modepy.spaces + +.. class:: FunctionSpace + + See :class:`modepy.FunctionSpace`. + +.. class:: PN + + See :class:`modepy.PN`. + +.. class:: QN + + See :class:`modepy.QN`. +""" + +__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 functools import singledispatch +from modepy.shapes import Shape, Simplex, Hypercube + + +# {{{ function spaces + +class FunctionSpace: + r"""An opaque object representing a finite-dimensional function space + of functions :math:`\mathbb R^n \to :math:`\mathbb R`. + + .. attribute:: spatial_dim + + :math:`n` in the above definition, the number of spatial dimensions + in which the functions in the space operate. + + .. attribute:: space_dim + + The number of dimensions of the function space. + """ + + +class PN(FunctionSpace): + r"""The function space of polynomials with total degree :math:`N`=:attr:`order`. + + .. math:: + + P^N:=\operatorname{span}\left\{\prod_{i=1}^d x_i^{n_i}:\sum n_i\le N\right\}. + + .. automethod:: __init__ + .. attribute:: order + """ + def __init__(self, spatial_dim, order): + super().__init__() + self.spatial_dim = spatial_dim + self.order = order + + @property + def space_dim(self): + spdim = self.spatial_dim + order = self.order + try: + from math import comb # comb is v3.8+ + return comb(order + spdim, spdim) + except ImportError: + from functools import reduce + from operator import mul + return reduce(mul, range(order + 1, order + spdim + 1), 1) \ + // reduce(mul, range(1, spdim + 1), 1) + + +class QN(FunctionSpace): + r"""The function space of polynomials with maximum degree + :math:`N`=:attr:`order`: + + .. math:: + + Q^N:=\operatorname{span} + \left \{\prod_{i=1}^d x_i^{n_i}:\max n_i\le N\right\}. + + .. automethod:: __init__ + .. attribute:: order + """ + def __init__(self, spatial_dim, order): + super().__init__() + self.spatial_dim = spatial_dim + self.order = order + + @property + def space_dim(self): + return (self.order + 1)**self.spatial_dim + + +@singledispatch +def space_for_shape(shape: Shape, order: int) -> FunctionSpace: + r"""Return an unspecified instance of :class:`FunctionSpace` suitable + for approximation on *shape* attaining interpolation error of + :math:`O(h^{\text{order}+1})`. + """ + raise NotImplementedError(type(shape).__name__) + + +@space_for_shape.register(Simplex) +def _(shape: Simplex, order: int): + return PN(shape.dim, order) + + +@space_for_shape.register(Hypercube) +def _(shape: Hypercube, order: int): + return QN(shape.dim, order) + +# }}} + + +# vim: foldmethod=marker diff --git a/modepy/tools.py b/modepy/tools.py index e41a9502..9c99da3b 100644 --- a/modepy/tools.py +++ b/modepy/tools.py @@ -1,3 +1,21 @@ +""" +Transformations between coordinate systems on the simplex +--------------------------------------------------------- + +All of these expect and return arrays of shape *(dims, npts)*. + +.. autofunction:: equilateral_to_unit +.. autofunction:: barycentric_to_unit +.. autofunction:: unit_to_barycentric +.. autofunction:: barycentric_to_equilateral + +Interpolation quality +--------------------- + +.. autofunction:: estimate_lebesgue_constant + +""" + __copyright__ = "Copyright (C) 2013 Andreas Kloeckner" __license__ = """ @@ -26,6 +44,7 @@ import numpy.linalg as la from math import sqrt from pytools import memoize_method, MovedFunctionDeprecationWrapper +import modepy.shapes as shp try: @@ -146,6 +165,10 @@ def inverse(self): """The inverse :class:`AffineMap` of *self*.""" return AffineMap(la.inv(self.a), -la.solve(self.a, self.b)) +# }}} + + +# {{{ simplex coordinate mapping EQUILATERAL_TO_UNIT_MAP = { 1: AffineMap([[1]], [0]), @@ -166,30 +189,21 @@ def equilateral_to_unit(equi): def unit_vertices(dim): - result = np.empty((dim+1, dim), np.float64) - result.fill(-1) + from warnings import warn + warn("unit_vertices is deprecated. " + "Use modepy.biunit_vertices_for_shape instead. " + "unit_vertices will go away in 2022.", + DeprecationWarning, stacklevel=2) - for i in range(dim): - result[i+1, i] = 1 - - return result - - -# this should go away -UNIT_VERTICES = { - 0: unit_vertices(0), - 1: unit_vertices(1), - 2: unit_vertices(2), - 3: unit_vertices(3), - } + return shp.biunit_vertices_for_shape(shp.Simplex(dim)).T def barycentric_to_unit(bary): """ - :arg bary: shaped ``(dims+1,npoints)`` + :arg bary: shaped ``(dims+1, npoints)`` """ dims = len(bary)-1 - return np.dot(unit_vertices(dims).T, bary) + return np.dot(shp.biunit_vertices_for_shape(shp.Simplex(dims)), bary) def unit_to_barycentric(unit): @@ -229,23 +243,7 @@ def barycentric_to_equilateral(bary): # }}} -def pick_random_simplex_unit_coordinate(rng, dims): - offset = 0.05 - base = -1 + offset - remaining = 2 - dims*offset - r = np.zeros(dims, np.float64) - for j in range(dims): - rn = rng.uniform(0, remaining) - r[j] = base + rn - remaining -= rn - return r - - -def pick_random_hypercube_unit_coordinate(rng, dims): - return np.array([rng.uniform(-1.0, 1.0) for _ in range(dims)]) - - -# {{{ submeshes, plotting helpers +# {{{ submeshes def simplex_submesh(node_tuples): """Return a list of tuples of indices into the node list that @@ -258,88 +256,7 @@ def simplex_submesh(node_tuples): :func:`pytools.generate_nonnegative_integer_tuples_summing_to_at_most` may be used to generate *node_tuples*. """ - from pytools import single_valued, add_tuples - dims = single_valued(len(nt) for nt in node_tuples) - - node_dict = { - ituple: idx - for idx, ituple in enumerate(node_tuples)} - - if dims == 1: - result = [] - - def try_add_line(d1, d2): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - )) - except KeyError: - pass - - for current in node_tuples: - try_add_line((0,), (1,),) - - return result - elif dims == 2: - # {{{ triangle sub-mesh - result = [] - - def try_add_tri(d1, d2, d3): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - )) - except KeyError: - pass - - for current in node_tuples: - # this is a tesselation of a square into two triangles. - # subtriangles that fall outside of the master triangle are - # simply not added. - - # positively oriented - try_add_tri((0, 0), (1, 0), (0, 1)) - try_add_tri((1, 0), (1, 1), (0, 1)) - - return result - - # }}} - elif dims == 3: - # {{{ tet sub-mesh - - def try_add_tet(d1, d2, d3, d4): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - node_dict[add_tuples(current, d4)], - )) - except KeyError: - pass - - result = [] - for current in node_tuples: - # this is a tesselation of a cube into six tets. - # subtets that fall outside of the master tet are simply not added. - - # positively oriented - try_add_tet((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) - try_add_tet((1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0)) - try_add_tet((1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1)) - - try_add_tet((1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0)) - try_add_tet((0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1)) - try_add_tet((0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0)) - - return result - - # }}} - else: - raise NotImplementedError("%d-dimensional sub-meshes" % dims) + return shp.submesh_for_shape(shp.Simplex(len(node_tuples[0])), node_tuples) submesh = MovedFunctionDeprecationWrapper(simplex_submesh) @@ -360,29 +277,18 @@ def hypercube_submesh(node_tuples): .. versionadded:: 2020.2 """ + from warnings import warn + warn("hypercube_submesh is deprecated. " + "Use submesh_for_shape instead. " + "hypercube_submesh will go away in 2022.", + DeprecationWarning, stacklevel=2) - from pytools import single_valued, add_tuples - dims = single_valued(len(nt) for nt in node_tuples) - - node_dict = { - ituple: idx - for idx, ituple in enumerate(node_tuples)} - - from pytools import generate_nonnegative_integer_tuples_below as gnitb - - result = [] - for current in node_tuples: - try: - result.append(tuple( - node_dict[add_tuples(current, offset)] - for offset in gnitb(2, dims))) - - except KeyError: - pass + return shp.submesh_for_shape(shp.Hypercube(len(node_tuples[0])), node_tuples) - return result +# }}} +# {{{ plotting helpers def plot_element_values(n, nodes, values, resample_n=None, node_tuples=None, show_nodes=False): dims = len(nodes) @@ -424,32 +330,23 @@ def plot_element_values(n, nodes, values, resample_n=None, # {{{ lebesgue constant -def _evaluate_lebesgue_function(n, nodes, domain): - dims = len(nodes) +def _evaluate_lebesgue_function(n, nodes, shape): huge_n = 30*n - if domain == "simplex": - from modepy.modes import simplex_onb as domain_basis_onb - from pytools import ( - generate_nonnegative_integer_tuples_summing_to_at_most - as generate_node_tuples) - elif domain == "hypercube": - from modepy.modes import ( - legendre_tensor_product_basis as domain_basis_onb) - from pytools import ( - generate_nonnegative_integer_tuples_below - as generate_node_tuples) - else: - raise ValueError(f"unknown domain: '{domain}'") + from modepy.spaces import space_for_shape + from modepy.modes import basis_for_space + from modepy.nodes import node_tuples_for_space + space = space_for_shape(shape, n) + huge_space = space_for_shape(shape, huge_n) - basis = domain_basis_onb(dims, n) - equi_node_tuples = list(generate_node_tuples(huge_n, dims)) + basis = basis_for_space(space) + equi_node_tuples = node_tuples_for_space(huge_space) equi_nodes = (np.array(equi_node_tuples, dtype=np.float64)/huge_n*2 - 1).T from modepy.matrices import vandermonde - vdm = vandermonde(basis, nodes) + vdm = vandermonde(basis.functions, nodes) - eq_vdm = vandermonde(basis, equi_nodes) + eq_vdm = vandermonde(basis.functions, equi_nodes) eq_to_out = la.solve(vdm.T, eq_vdm.T).T lebesgue_worst = np.sum(np.abs(eq_to_out), axis=1) @@ -457,16 +354,15 @@ def _evaluate_lebesgue_function(n, nodes, domain): return lebesgue_worst, equi_node_tuples, equi_nodes -def estimate_lebesgue_constant(n, nodes, domain=None, visualize=False): +def estimate_lebesgue_constant(n, nodes, shape=None, visualize=False): """Estimate the `Lebesgue constant `_ of the *nodes* at polynomial order *n*. - :arg nodes: an array of shape *(dims, nnodes)* as returned by + :arg nodes: an array of shape *(dim, nnodes)* as returned by :func:`modepy.warp_and_blend_nodes`. - :arg domain: represents the domain of the reference element and can be - either ``"simplex"`` or ``"hypercube"``. + :arg shape: a :class:`~modepy.shapes.Shape`. :arg visualize: visualize the function that gives rise to the returned Lebesgue constant. (2D only for now) :return: the Lebesgue constant, a scalar. @@ -477,27 +373,32 @@ def estimate_lebesgue_constant(n, nodes, domain=None, visualize=False): *domain* parameter was added with support for nodes on the unit hypercube (i.e. unit square in 2D and unit cube in 3D). + + .. versionchanged:: 2020.3 + + Renamed *domain* to *shape*. """ - if domain is None: - domain = "simplex" + dim = len(nodes) + if shape is None: + from warnings import warn + warn("Not passing shape is deprecated and will stop working " + "in 2022.", DeprecationWarning, stacklevel=2) + from modepy.shapes import Simplex + shape = Simplex(dim) + else: + if shape.dim != dim: + raise ValueError(f"expected {shape.dim}-dimensional nodes") - dims = len(nodes) lebesgue_worst, equi_node_tuples, equi_nodes = \ - _evaluate_lebesgue_function(n, nodes, domain) + _evaluate_lebesgue_function(n, nodes, shape) lebesgue_constant = np.max(lebesgue_worst) if not visualize: return lebesgue_constant - if dims == 2: + if shape.dim == 2: print(f"Lebesgue constant: {lebesgue_constant}") - - if domain == "simplex": - triangles = simplex_submesh(equi_node_tuples) - elif domain == "hypercube": - triangles = hypercube_submesh(equi_node_tuples) - else: - triangles = None + triangles = shp.submesh_for_shape(shape, equi_node_tuples) try: import mayavi.mlab as mlab @@ -532,7 +433,7 @@ def estimate_lebesgue_constant(n, nodes, domain=None, visualize=False): ax.set_aspect("equal") plt.show() else: - raise ValueError(f"visualization is not supported in {dims}D") + raise ValueError(f"visualization is not supported in {shape.dim}D") return lebesgue_constant diff --git a/modepy/version.py b/modepy/version.py index 82ea6fb1..acb402a2 100644 --- a/modepy/version.py +++ b/modepy/version.py @@ -1,2 +1,2 @@ -VERSION = (2020, 2) +VERSION = (2020, 3) VERSION_TEXT = ".".join(str(i) for i in VERSION) diff --git a/setup.py b/setup.py index 337a13a8..05c8d197 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,8 @@ def main(): "numpy", "pytools>=2013.1", "pytest>=2.3", + + "dataclasses; python_version<='3.6'", ]) diff --git a/test/test_modes.py b/test/test_modes.py index ffe61906..f50ca904 100644 --- a/test/test_modes.py +++ b/test/test_modes.py @@ -21,10 +21,11 @@ """ +from functools import partial import numpy as np import numpy.linalg as la import pytest -import modepy.modes as m +import modepy as mp from pymbolic.mapper.stringifier import ( CSESplittingStringifyMapperMixin, StringifyMapper) from pymbolic.mapper.evaluator import EvaluationMapper @@ -51,7 +52,7 @@ def test_orthonormality_jacobi_1d(alpha, beta, ebound): quad = JacobiGaussQuadrature(alpha, beta, 4*max_n) from functools import partial - jac_f = [partial(m.jacobi, alpha, beta, n) for n in range(max_n)] + jac_f = [partial(mp.jacobi, alpha, beta, n) for n in range(max_n)] maxerr = 0 for i, fi in enumerate(jac_f): @@ -76,19 +77,22 @@ def test_orthonormality_jacobi_1d(alpha, beta, ebound): # (7, 3e-14), # (9, 2e-13), ]) -@pytest.mark.parametrize("dims", [2, 3]) -def test_pkdo_orthogonality(dims, order, ebound): - """Test orthogonality of simplicial bases using Grundmann-Moeller cubature.""" - - from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature - from modepy.modes import simplex_onb +@pytest.mark.parametrize("shape", [ + mp.Simplex(2), + mp.Simplex(3), + mp.Hypercube(2), + mp.Hypercube(3), + ]) +def test_orthogonality(shape, order, ebound): + """Test orthogonality of ONBs using cubature.""" - cub = GrundmannMoellerSimplexQuadrature(order, dims) - basis = simplex_onb(dims, order) + qspace = mp.space_for_shape(shape, 2*order) + cub = mp.quadrature_for_space(qspace, shape) + basis = mp.orthonormal_basis_for_space(mp.space_for_shape(shape, order)) maxerr = 0 - for i, f in enumerate(basis): - for j, g in enumerate(basis): + for i, f in enumerate(basis.functions): + for j, g in enumerate(basis.functions): if i == j: true_result = 1 else: @@ -103,52 +107,70 @@ def test_pkdo_orthogonality(dims, order, ebound): # print(order, maxerr) -@pytest.mark.parametrize("dims", [1, 2, 3]) +def get_inhomogeneous_tensor_prod_basis(space): + # FIXME: Yuck. A total lie. Not a basis for the space at all. + assert isinstance(space, mp.QN) + orders = (3, 5, 7)[:space.spatial_dim] + + return mp.TensorProductBasis( + [[partial(mp.jacobi, 0, 0, n) for n in range(o)] + for o in orders], + [[partial(mp.grad_jacobi, 0, 0, n) for n in range(o)] + for o in orders], + orth_weight=1) + + +@pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("order", [5, 8]) -@pytest.mark.parametrize(("eltype", "basis_getter", "grad_basis_getter"), [ - ("simplex", m.simplex_onb, m.grad_simplex_onb), - ("simplex", m.simplex_monomial_basis, m.grad_simplex_monomial_basis), - ("tensor", m.legendre_tensor_product_basis, m.grad_legendre_tensor_product_basis) +@pytest.mark.parametrize(("shape_cls", "basis_getter"), [ + (mp.Simplex, mp.basis_for_space), + (mp.Simplex, mp.orthonormal_basis_for_space), + (mp.Simplex, mp.monomial_basis_for_space), + + (mp.Hypercube, mp.basis_for_space), + (mp.Hypercube, mp.orthonormal_basis_for_space), + (mp.Hypercube, mp.monomial_basis_for_space), + + (mp.Hypercube, get_inhomogeneous_tensor_prod_basis), ]) -def test_basis_grad(dims, order, eltype, basis_getter, grad_basis_getter): +def test_basis_grad(dim, shape_cls, order, basis_getter): """Do a simplistic FD-style check on the gradients of the basis.""" h = 1.0e-4 - if eltype == "simplex" and order == 8 and dims == 3: - factor = 3.0 - else: - factor = 1.0 - - if eltype == "simplex": - from modepy.tools import \ - pick_random_simplex_unit_coordinate as pick_random_unit_coordinate - elif eltype == "tensor": - from modepy.tools import \ - pick_random_hypercube_unit_coordinate as pick_random_unit_coordinate - else: - raise ValueError(f"unknown element type: {eltype}") - from random import Random - rng = Random(17) + shape = shape_cls(dim) + rng = np.random.Generator(np.random.PCG64(17)) + basis = basis_getter(mp.space_for_shape(shape, order)) + from pytools.convergence import EOCRecorder from pytools import wandering_element for i_bf, (bf, gradbf) in enumerate(zip( - basis_getter(dims, order), - grad_basis_getter(dims, order), + basis.functions, + basis.gradients, )): - for i in range(10): - r = pick_random_unit_coordinate(rng, dims) + eoc_rec = EOCRecorder() + for h in [1e-2, 1e-3]: + r = mp.random_nodes_for_shape(shape, nnodes=1000, rng=rng) gradbf_v = np.array(gradbf(r)) gradbf_v_num = np.array([ (bf(r+h*unit) - bf(r-h*unit))/(2*h) - for unit_tuple in wandering_element(dims) - for unit in (np.array(unit_tuple),) + for unit_tuple in wandering_element(shape.dim) + for unit in (np.array(unit_tuple).reshape(-1, 1),) ]) - err = la.norm(gradbf_v_num - gradbf_v) + ref_norm = la.norm((gradbf_v).reshape(-1), np.inf) + err = la.norm((gradbf_v_num - gradbf_v).reshape(-1), np.inf) + if ref_norm > 1e-13: + err = err/ref_norm + logger.info("error: %.5", err) - assert err < factor * h, (err, i_bf) + eoc_rec.add_data_point(h, err) + + tol = 1e-8 + if eoc_rec.max_error() >= tol: + print(eoc_rec) + assert (eoc_rec.max_error() < tol or eoc_rec.order_estimate() >= 1.5) # {{{ test symbolic modes @@ -163,20 +185,23 @@ def map_if(self, expr): self.rec(expr.then), self.rec(expr.else_)) -@pytest.mark.parametrize(("domain", "get_basis", "get_grad_basis"), [ - ("simplex", m.simplex_onb, m.grad_simplex_onb), - ("simplex", m.simplex_monomial_basis, m.grad_simplex_monomial_basis), - ("hypercube", m.legendre_tensor_product_basis, - m.grad_legendre_tensor_product_basis), +@pytest.mark.parametrize("shape", [ + mp.Simplex(1), + mp.Simplex(2), + mp.Simplex(3), + mp.Hypercube(1), + mp.Hypercube(2), + mp.Hypercube(3), ]) -@pytest.mark.parametrize("dims", [1, 2, 3]) -@pytest.mark.parametrize("n", [5, 10]) -def test_symbolic_basis(domain, dims, n, - get_basis, - get_grad_basis): - - basis = get_basis(dims, n) - sym_basis = [m.symbolicize_function(f, dims) for f in basis] +@pytest.mark.parametrize("order", [5, 8]) +@pytest.mark.parametrize("basis_getter", [ + (mp.basis_for_space), + (mp.orthonormal_basis_for_space), + (mp.monomial_basis_for_space), + ]) +def test_symbolic_basis(shape, order, basis_getter): + basis = basis_getter(mp.space_for_shape(shape, order)) + sym_basis = [mp.symbolicize_function(f, shape.dim) for f in basis.functions] # {{{ test symbolic against direct eval @@ -184,15 +209,10 @@ def test_symbolic_basis(domain, dims, n, print("VALUES") print(75*"#") - r = np.random.rand(dims, 10000)*2 - 1 - if domain == "simplex": - r = r[:, r.sum(axis=0) < 0] - elif domain == "hypercube": - pass - else: - raise ValueError(f"unexpected domain: {domain}") + rng = np.random.Generator(np.random.PCG64(17)) + r = mp.random_nodes_for_shape(shape, 10000, rng=rng) - for func, sym_func in zip(basis, sym_basis): + for func, sym_func in zip(basis.functions, sym_basis): strmap = MyStringifyMapper() s = strmap(sym_func) for name, val in strmap.cse_name_list: @@ -218,10 +238,9 @@ def test_symbolic_basis(domain, dims, n, print("GRADIENTS") print(75*"#") - grad_basis = get_grad_basis(dims, n) - sym_grad_basis = [m.symbolicize_function(f, dims) for f in grad_basis] + sym_grad_basis = [mp.symbolicize_function(f, shape.dim) for f in basis.gradients] - for grad, sym_grad in zip(grad_basis, sym_grad_basis): + for grad, sym_grad in zip(basis.gradients, sym_grad_basis): strmap = MyStringifyMapper() s = strmap(sym_grad) for name, val in strmap.cse_name_list: diff --git a/test/test_nodes.py b/test/test_nodes.py index c7b52092..05f6941c 100644 --- a/test/test_nodes.py +++ b/test/test_nodes.py @@ -24,24 +24,21 @@ import numpy as np import numpy.linalg as la import pytest +import modepy.shapes as shp +import modepy.nodes as nd @pytest.mark.parametrize("dims", [1, 2, 3]) def test_barycentric_coordinate_map(dims): - from random import Random - rng = Random(17) - - n = 5 - unit = np.empty((dims, n)) + n = 100 from modepy.tools import ( - pick_random_simplex_unit_coordinate, unit_to_barycentric, barycentric_to_unit, barycentric_to_equilateral, equilateral_to_unit,) - for i in range(n): - unit[:, i] = pick_random_simplex_unit_coordinate(rng, dims) + rng = np.random.Generator(np.random.PCG64(17)) + unit = nd.random_nodes_for_shape(shp.Simplex(dims), n, rng=rng) bary = unit_to_barycentric(unit) assert (np.abs(np.sum(bary, axis=0) - 1) < 1e-15).all() @@ -57,11 +54,10 @@ def test_barycentric_coordinate_map(dims): def test_warp(): """Check some assumptions on the node warp factor calculator""" n = 17 - from modepy.nodes import warp_factor from functools import partial def wfc(x): - return warp_factor(n, np.array([x]), scaled=False)[0] + return nd.warp_factor(n, np.array([x]), scaled=False)[0] assert abs(wfc(-1)) < 1e-12 assert abs(wfc(1)) < 1e-12 @@ -69,7 +65,7 @@ def wfc(x): from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature lgq = LegendreGaussQuadrature(n) - assert abs(lgq(partial(warp_factor, n, scaled=False))) < 6e-14 + assert abs(lgq(partial(nd.warp_factor, n, scaled=False))) < 6e-14 def test_tri_face_node_distribution(): @@ -85,8 +81,7 @@ def test_tri_face_node_distribution(): as gnitstam node_tuples = list(gnitstam(n, 2)) - from modepy.nodes import warp_and_blend_nodes - unodes = warp_and_blend_nodes(2, n, node_tuples) + unodes = nd.warp_and_blend_nodes(2, n, node_tuples) faces = [ [i for i, nt in enumerate(node_tuples) if nt[0] == 0], @@ -116,8 +111,7 @@ def test_simp_nodes(dims, n): eps = 1e-10 - from modepy.nodes import warp_and_blend_nodes - unodes = warp_and_blend_nodes(dims, n) + unodes = nd.warp_and_blend_nodes(dims, n) assert (unodes >= -1-eps).all() assert (np.sum(unodes) <= eps).all() @@ -139,12 +133,11 @@ def test_affine_map(): @pytest.mark.parametrize("dim", [1, 2, 3, 4]) def test_tensor_product_nodes(dim): - from modepy.nodes import tensor_product_nodes nnodes = 10 nodes_1d = np.arange(nnodes) - nodes = tensor_product_nodes(dim, nodes_1d) + nodes = nd.tensor_product_nodes(dim, nodes_1d) assert np.allclose( - nodes[-1], + nodes[0], np.array(nodes_1d.tolist() * nnodes**(dim - 1))) diff --git a/test/test_quadrature.py b/test/test_quadrature.py index dcfa0e27..8a3912df 100644 --- a/test/test_quadrature.py +++ b/test/test_quadrature.py @@ -23,6 +23,7 @@ import numpy as np import pytest + import modepy as mp import logging @@ -144,11 +145,11 @@ def test_simplex_quadrature(quad_class, highest_order, dim): break -@pytest.mark.parametrize("cls", [ +@pytest.mark.parametrize("quad_class", [ mp.WitherdenVincentQuadrature ]) @pytest.mark.parametrize("dim", [2, 3]) -def test_hypercube_quadrature(cls, dim): +def test_hypercube_quadrature(quad_class, dim): from pytools import \ generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam from modepy.tools import Monomial @@ -167,14 +168,15 @@ def _check_monomial(quad, comb): order = 1 while True: try: - quad = cls(order, dim) + quad = quad_class(order, dim) except mp.QuadratureRuleUnavailable: logger.info("UNAVAILABLE at order %d", order) break assert np.all(quad.weights > 0) - logger.info("quadrature: %s %d %d", cls, order, quad.exact_to) + logger.info("quadrature: %s %d %d", + quad_class.__name__.lower(), order, quad.exact_to) for comb in gnitstam(quad.exact_to, dim): assert _check_monomial(quad, comb) < 5.0e-15 diff --git a/test/test_tools.py b/test/test_tools.py index b65e7b61..b1921171 100644 --- a/test/test_tools.py +++ b/test/test_tools.py @@ -88,9 +88,10 @@ def constant(x): ("c1-2d", partial(c1, -0.1), 2, 15, -2.3), ]) def test_modal_decay(case_name, test_func, dims, n, expected_expn): + space = mp.PN(dims, n) nodes = mp.warp_and_blend_nodes(dims, n) - basis = mp.simplex_onb(dims, n) - vdm = mp.vandermonde(basis, nodes) + basis = mp.orthonormal_basis_for_space(space) + vdm = mp.vandermonde(basis.functions, nodes) f = test_func(nodes[0]) coeffs = la.solve(vdm, f) @@ -120,8 +121,8 @@ def test_modal_decay(case_name, test_func, dims, n, expected_expn): def test_residual_estimation(case_name, test_func, dims, n): def estimate_resid(inner_n): nodes = mp.warp_and_blend_nodes(dims, inner_n) - basis = mp.simplex_onb(dims, inner_n) - vdm = mp.vandermonde(basis, nodes) + basis = mp.orthonormal_basis_for_space(mp.PN(dims, inner_n)) + vdm = mp.vandermonde(basis.functions, nodes) f = test_func(nodes[0]) coeffs = la.solve(vdm, f) @@ -141,36 +142,29 @@ def estimate_resid(inner_n): # {{{ test_resampling_matrix @pytest.mark.parametrize("dims", [1, 2, 3]) -@pytest.mark.parametrize("eltype", ["simplex", "tensor"]) -def test_resampling_matrix(dims, eltype): - ncoarse = 5 - nfine = 10 - - if eltype == "simplex": - coarse_nodes = mp.warp_and_blend_nodes(dims, ncoarse) - fine_nodes = mp.warp_and_blend_nodes(dims, nfine) - - coarse_basis = mp.simplex_onb(dims, ncoarse) - fine_basis = mp.simplex_onb(dims, nfine) - elif eltype == "tensor": - coarse_nodes = mp.legendre_gauss_lobatto_tensor_product_nodes(dims, ncoarse) - fine_nodes = mp.legendre_gauss_lobatto_tensor_product_nodes(dims, nfine) - - coarse_basis = mp.legendre_tensor_product_basis(dims, ncoarse) - fine_basis = mp.legendre_tensor_product_basis(dims, nfine) - else: - raise ValueError(f"unknown element type: {eltype}") +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) +def test_resampling_matrix(dims, shape_cls, ncoarse=5, nfine=10): + shape = shape_cls(dims) + + coarse_space = mp.space_for_shape(shape, ncoarse) + fine_space = mp.space_for_shape(shape, nfine) + + coarse_nodes = mp.edge_clustered_nodes_for_space(coarse_space, shape) + coarse_basis = mp.basis_for_space(coarse_space) + + fine_nodes = mp.edge_clustered_nodes_for_space(fine_space, shape) + fine_basis = mp.basis_for_space(fine_space) my_eye = np.dot( - mp.resampling_matrix(fine_basis, coarse_nodes, fine_nodes), - mp.resampling_matrix(coarse_basis, fine_nodes, coarse_nodes)) + mp.resampling_matrix(fine_basis.functions, coarse_nodes, fine_nodes), + mp.resampling_matrix(coarse_basis.functions, fine_nodes, coarse_nodes)) assert la.norm(my_eye - np.eye(len(my_eye))) < 3e-13 my_eye_least_squares = np.dot( - mp.resampling_matrix(coarse_basis, coarse_nodes, fine_nodes, + mp.resampling_matrix(coarse_basis.functions, coarse_nodes, fine_nodes, least_squares_ok=True), - mp.resampling_matrix(coarse_basis, fine_nodes, coarse_nodes), + mp.resampling_matrix(coarse_basis.functions, fine_nodes, coarse_nodes), ) assert la.norm(my_eye_least_squares - np.eye(len(my_eye_least_squares))) < 4e-13 @@ -181,22 +175,14 @@ def test_resampling_matrix(dims, eltype): # {{{ test_diff_matrix @pytest.mark.parametrize("dims", [1, 2, 3]) -@pytest.mark.parametrize("eltype", ["simplex", "tensor"]) -def test_diff_matrix(dims, eltype): - n = 5 - - if eltype == "simplex": - nodes = mp.warp_and_blend_nodes(dims, n) - basis = mp.simplex_onb(dims, n) - grad_basis = mp.grad_simplex_onb(dims, n) - elif eltype == "tensor": - nodes = mp.legendre_gauss_lobatto_tensor_product_nodes(dims, n) - basis = mp.legendre_tensor_product_basis(dims, n) - grad_basis = mp.grad_legendre_tensor_product_basis(dims, n) - else: - raise ValueError(f"unknown element type: {eltype}") - - diff_mat = mp.differentiation_matrices(basis, grad_basis, nodes) +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) +def test_diff_matrix(dims, shape_cls, order=5): + shape = shape_cls(dims) + space = mp.space_for_shape(shape, order) + nodes = mp.edge_clustered_nodes_for_space(space, shape) + basis = mp.basis_for_space(space) + + diff_mat = mp.differentiation_matrices(basis.functions, basis.gradients, nodes) if isinstance(diff_mat, tuple): diff_mat = diff_mat[0] @@ -213,15 +199,16 @@ def test_diff_matrix(dims, eltype): @pytest.mark.parametrize("dims", [2, 3]) def test_diff_matrix_permutation(dims): order = 5 + space = mp.PN(dims, order) from pytools import \ generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam node_tuples = list(gnitstam(order, dims)) - simplex_onb = mp.simplex_onb(dims, order) - grad_simplex_onb = mp.grad_simplex_onb(dims, order) + simplex_onb = mp.orthonormal_basis_for_space(space) nodes = np.array(mp.warp_and_blend_nodes(dims, order, node_tuples=node_tuples)) - diff_matrices = mp.differentiation_matrices(simplex_onb, grad_simplex_onb, nodes) + diff_matrices = mp.differentiation_matrices( + simplex_onb.functions, simplex_onb.gradients, nodes) for iref_axis in range(dims): perm = mp.diff_matrix_permutation(node_tuples, iref_axis) @@ -233,65 +220,141 @@ def test_diff_matrix_permutation(dims): # }}} -# {{{ test_face_mass_matrix - -@pytest.mark.parametrize("dim", [1, 2, 3]) -def test_modal_face_mass_matrix(dim, order=3): - from modepy.tools import unit_vertices - all_verts = unit_vertices(dim).T +# {{{ face mass matrices (deprecated) - basis = mp.simplex_onb(dim, order) +@pytest.mark.parametrize("dims", [2, 3]) +def test_deprecated_modal_face_mass_matrix(dims, order=3): + # FIXME DEPRECATED remove along with modal_face_mass_matrix (>=2022) + shape = mp.Simplex(dims) + space = mp.space_for_shape(shape, order) - # np.set_printoptions(linewidth=200) + vertices = mp.biunit_vertices_for_shape(shape) + basis = mp.basis_for_space(space) from modepy.matrices import modal_face_mass_matrix - for iface in range(dim+1): - verts = np.hstack([all_verts[:, :iface], all_verts[:, iface+1:]]) + for face in mp.faces_for_shape(shape): + face_vertices = vertices[:, face.volume_vertex_indices] - fmm = modal_face_mass_matrix(basis, order, verts) - fmm2 = modal_face_mass_matrix(basis, order+1, verts) + fmm = modal_face_mass_matrix( + basis.functions, order, face_vertices) + fmm2 = modal_face_mass_matrix( + basis.functions, order+1, face_vertices) - assert la.norm(fmm-fmm2, np.inf) < 1e-11 + error = la.norm(fmm - fmm2, np.inf) / la.norm(fmm2, np.inf) + logger.info("fmm error: %.5e", error) + assert error < 1e-11, f"error {error:.5e} on face {face.face_index}" fmm[np.abs(fmm) < 1e-13] = 0 - - print(fmm) nnz = np.sum(fmm > 0) - print(nnz) + logger.info("fmm: nnz %d\n%s", nnz, fmm) -@pytest.mark.parametrize("dim", [1, 2, 3]) -def test_nodal_face_mass_matrix(dim, order=3): - from modepy.tools import unit_vertices - all_verts = unit_vertices(dim).T - basis = mp.simplex_onb(dim, order) +@pytest.mark.parametrize("dims", [2, 3]) +def test_deprecated_nodal_face_mass_matrix(dims, order=3): + # FIXME DEPRECATED remove along with nodal_face_mass_matrix (>=2022) + vol_shape = mp.Simplex(dims) + vol_space = mp.space_for_shape(vol_shape, order) - np.set_printoptions(linewidth=200) + vertices = mp.biunit_vertices_for_shape(vol_shape) + volume_nodes = mp.edge_clustered_nodes_for_space(vol_space, vol_shape) + volume_basis = mp.basis_for_space(vol_space) from modepy.matrices import nodal_face_mass_matrix - volume_nodes = mp.warp_and_blend_nodes(dim, order) - face_nodes = mp.warp_and_blend_nodes(dim-1, order) + for face in mp.faces_for_shape(vol_shape): + face_space = mp.space_for_shape(face, order) + face_nodes = mp.edge_clustered_nodes_for_space(face_space, face) + face_vertices = vertices[:, face.volume_vertex_indices] + + fmm = nodal_face_mass_matrix( + volume_basis.functions, volume_nodes, + face_nodes, order, face_vertices) + fmm2 = nodal_face_mass_matrix( + volume_basis.functions, + volume_nodes, face_nodes, order+1, face_vertices) + + error = la.norm(fmm - fmm2, np.inf) / la.norm(fmm2, np.inf) + logger.info("fmm error: %.5e", error) + assert error < 5e-11, f"error {error:.5e} on face {face.face_index}" - for iface in range(dim+1): - verts = np.hstack([all_verts[:, :iface], all_verts[:, iface+1:]]) + fmm[np.abs(fmm) < 1e-13] = 0 + nnz = np.sum(fmm > 0) + + logger.info("fmm: nnz %d\n%s", nnz, fmm) - fmm = nodal_face_mass_matrix(basis, volume_nodes, face_nodes, order, - verts) - fmm2 = nodal_face_mass_matrix(basis, volume_nodes, face_nodes, order+1, - verts) + logger.info("mass matrix:\n%s", mp.mass_matrix( + mp.basis_for_space(face_space).functions, + mp.edge_clustered_nodes_for_space(face_space, face))) - assert la.norm(fmm-fmm2, np.inf) < 1e-11 +# }}} + + +# {{{ face mass matrices + +@pytest.mark.parametrize("dims", [2, 3]) +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) +def test_modal_mass_matrix_for_face(dims, shape_cls, order=3): + vol_shape = shape_cls(dims) + vol_space = mp.space_for_shape(vol_shape, order) + vol_basis = mp.basis_for_space(vol_space) + + from modepy.matrices import modal_mass_matrix_for_face + for face in mp.faces_for_shape(vol_shape): + face_space = mp.space_for_shape(face, order) + face_basis = mp.basis_for_space(face_space) + face_quad = mp.quadrature_for_space(mp.space_for_shape(face, 2*order), face) + face_quad2 = mp.quadrature_for_space( + mp.space_for_shape(face, 2*order+2), face) + fmm = modal_mass_matrix_for_face( + face, face_quad, face_basis.functions, vol_basis.functions) + fmm2 = modal_mass_matrix_for_face( + face, face_quad2, face_basis.functions, vol_basis.functions) + + error = la.norm(fmm - fmm2, np.inf) / la.norm(fmm2, np.inf) + logger.info("fmm error: %.5e", error) + assert error < 1e-11, f"error {error:.5e} on face {face.face_index}" fmm[np.abs(fmm) < 1e-13] = 0 + nnz = np.sum(fmm > 0) + + logger.info("fmm: nnz %d\n%s", nnz, fmm) - print(fmm) - nnz = np.sum(np.abs(fmm) > 0) - print(nnz) - print(mp.mass_matrix( - mp.simplex_onb(dim-1, order), - mp.warp_and_blend_nodes(dim-1, order), )) +@pytest.mark.parametrize("dims", [2, 3]) +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) +def test_nodal_mass_matrix_for_face(dims, shape_cls, order=3): + vol_shape = shape_cls(dims) + vol_space = mp.space_for_shape(vol_shape, order) + + volume_nodes = mp.edge_clustered_nodes_for_space(vol_space, vol_shape) + volume_basis = mp.basis_for_space(vol_space) + + from modepy.matrices import nodal_mass_matrix_for_face + for face in mp.faces_for_shape(vol_shape): + face_space = mp.space_for_shape(face, order) + face_basis = mp.basis_for_space(face_space) + face_nodes = mp.edge_clustered_nodes_for_space(face_space, face) + face_quad = mp.quadrature_for_space(mp.space_for_shape(face, 2*order), face) + face_quad2 = mp.quadrature_for_space( + mp.space_for_shape(face, 2*order+2), face) + fmm = nodal_mass_matrix_for_face( + face, face_quad, face_basis.functions, volume_basis.functions, + volume_nodes, face_nodes) + fmm2 = nodal_mass_matrix_for_face( + face, face_quad2, face_basis.functions, volume_basis.functions, + volume_nodes, face_nodes) + + error = la.norm(fmm - fmm2, np.inf) / la.norm(fmm2, np.inf) + logger.info("fmm error: %.5e", error) + assert error < 5e-11, f"error {error:.5e} on face {face.face_index}" + + fmm[np.abs(fmm) < 1e-13] = 0 + nnz = np.sum(fmm > 0) + + logger.info("fmm: nnz %d\n%s", nnz, fmm) + + logger.info("mass matrix:\n%s", + mp.mass_matrix(face_basis.functions, face_nodes)) # }}} @@ -300,28 +363,24 @@ def test_nodal_face_mass_matrix(dim, order=3): @pytest.mark.parametrize("dims", [1, 2]) @pytest.mark.parametrize("order", [3, 5, 8]) -@pytest.mark.parametrize("domain", ["simplex", "hypercube"]) -def test_estimate_lebesgue_constant(dims, order, domain, visualize=False): +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) +def test_estimate_lebesgue_constant(dims, order, shape_cls, visualize=False): logging.basicConfig(level=logging.INFO) + shape = shape_cls(dims) + space = mp.space_for_shape(shape, order) - if domain == "simplex": - nodes = mp.warp_and_blend_nodes(dims, order) - elif domain == "hypercube": - from modepy.nodes import legendre_gauss_lobatto_tensor_product_nodes - nodes = legendre_gauss_lobatto_tensor_product_nodes(dims, order) - else: - raise ValueError(f"unknown domain: '{domain}'") + nodes = mp.edge_clustered_nodes_for_space(space, shape) from modepy.tools import estimate_lebesgue_constant - lebesgue_constant = estimate_lebesgue_constant(order, nodes, domain=domain) - logger.info("%s-%d/%s: %.5e", domain, dims, order, lebesgue_constant) + lebesgue_constant = estimate_lebesgue_constant(order, nodes, shape=shape) + logger.info("%s-%d/%s: %.5e", shape, dims, order, lebesgue_constant) if not visualize: return from modepy.tools import _evaluate_lebesgue_function lebesgue, equi_node_tuples, equi_nodes = \ - _evaluate_lebesgue_function(order, nodes, domain) + _evaluate_lebesgue_function(order, nodes, shape) import matplotlib.pyplot as plt fig = plt.figure() @@ -342,7 +401,8 @@ def test_estimate_lebesgue_constant(dims, order, domain, visualize=False): else: raise ValueError(f"unsupported dimension: {dims}") - fig.savefig(f"estimate_lebesgue_constant_{domain}_{dims}_order_{order}") + shape_name = shape_cls.__name__.lower() + fig.savefig(f"estimate_lebesgue_constant_{shape_name}_{dims}_order_{order}") # }}} @@ -351,8 +411,8 @@ def test_estimate_lebesgue_constant(dims, order, domain, visualize=False): @pytest.mark.parametrize("dims", [2, 3, 4]) def test_hypercube_submesh(dims): - from modepy.tools import hypercube_submesh from pytools import generate_nonnegative_integer_tuples_below as gnitb + shape = mp.Hypercube(dims) node_tuples = list(gnitb(3, dims)) @@ -361,7 +421,7 @@ def test_hypercube_submesh(dims): assert len(node_tuples) == 3**dims - elements = hypercube_submesh(node_tuples) + elements = mp.submesh_for_shape(shape, node_tuples) for e in elements: logger.info("element: %s", e)