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 daf21889..c6d60997 100644 --- a/doc/nodes.rst +++ b/doc/nodes.rst @@ -1,39 +1,4 @@ Interpolation Nodes =================== -Simplices -^^^^^^^^^ - -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 - -Also see :class:`modepy.VioreanuRokhlinSimplexQuadrature` if nodes on the -boundary are not required. - -Hypercubes -^^^^^^^^^^ - -Node sets for interpolation ---------------------------- - -.. currentmodule:: modepy - -.. autofunction:: tensor_product_nodes -.. autofunction:: legendre_gauss_lobatto_tensor_product_nodes +.. automodule:: modepy.nodes diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 6d9045f2..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 ---------------------------------------- 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 b951e888..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,21 +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) + 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) + 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 @@ -58,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", @@ -65,16 +90,27 @@ "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", @@ -83,7 +119,6 @@ "ClenshawCurtisQuadrature", "FejerQuadrature", "WitherdenVincentQuadrature", - "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature", ] from pytools import MovedFunctionDeprecationWrapper diff --git a/modepy/matrices.py b/modepy/matrices.py index 91a40b8d..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,34 +245,62 @@ def mass_matrix(basis, nodes): return la.inv(inverse_mass_matrix(basis, nodes)) -def modal_face_mass_matrix(trial_basis, order, face_vertices, - test_basis=None, shape=None): +def modal_mass_matrix_for_face(face: Face, face_quad: Quadrature, + trial_functions, test_functions): + """ + .. versionadded:: 2020.3 """ - :arg face_vertices: an array of shape ``(dims, nvertices)``. - :arg shape: a :class:`~modepy.shapes.Shape` that identifies the - reference face element. - .. versionadded :: 2016.1 + mapped_nodes = face.map_to_volume(face_quad.nodes) + + result = np.empty((len(test_functions), len(trial_functions))) + + 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)) + - .. versionchanged:: 2020.3 +# {{{ deprecated junk - Added *shape* parameter and support for :math:`[-1, 1]^d` domains. +def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None): + """ + :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 - if shape is None: - from modepy.shapes import Simplex - shape = Simplex(face_vertices.shape[0]) + 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.shapes import get_face_map, get_quadrature - face = type(shape)(shape.dims - 1) - fmap = get_face_map(shape, face_vertices) - quad = get_quadrature(face, order) + 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) @@ -278,43 +309,36 @@ def modal_face_mass_matrix(trial_basis, order, face_vertices, 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 def nodal_face_mass_matrix(trial_basis, volume_nodes, face_nodes, order, - face_vertices, test_basis=None, shape=None): + face_vertices, test_basis=None): """ :arg face_vertices: an array of shape ``(dims, nvertices)``. - :arg shape: a :class:`~modepy.shapes.Shape` that identifies the - reference face element. .. versionadded :: 2016.1 - - .. versionchanged:: 2020.3 - - Added *shape* parameter and support for :math:`[-1, 1]^d` domains. """ + 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 - if shape is None: - from modepy.shapes import Simplex - shape = Simplex(face_vertices.shape[0]) - - from modepy.shapes import get_face_map - fmap = get_face_map(shape, 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, shape=shape) + 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 8e7c7240..52680048 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -22,16 +22,30 @@ """ +from warnings import warn import sys from math import sqrt +from functools import singledispatch, partial + import numpy as np -from modepy.shapes import Simplex, Hypercube -from modepy.shapes import get_basis, get_grad_basis, get_basis_with_mode_ids +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 -__doc__ = """:mod:`modepy.modes` provides orthonormal bases and their -derivatives on unit simplices. +Basis Retrieval +--------------- + +.. autoexception:: BasisNotOrthonormal +.. autoclass:: Basis + +.. autofunction:: basis_for_space +.. autofunction:: orthonormal_basis_for_space +.. autofunction:: monomial_basis_for_space Jacobi polynomials ------------------ @@ -41,8 +55,19 @@ .. 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 @@ -57,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 @@ -86,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`. """ @@ -442,14 +455,7 @@ def diff_monomial(r, o): # }}} -# {{{ dimension-independent interface for simplices - -def zerod_basis(x): - if len(x.shape) == 1: - return 1 - else: - return np.ones(x.shape[1]) - +# {{{ 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 @@ -471,25 +477,17 @@ def simplex_onb_with_mode_ids(dims, n): .. versionadded:: 2018.1 """ - from modepy.shapes import get_node_tuples - shape = Simplex(dims) - - from functools import partial - if dims == 0: - mode_ids = get_node_tuples(shape, n) - return mode_ids, (zerod_basis,) - elif dims == 1: - # FIXME: should also use get_node_tuples + warn("simplex_onb_with_mode_ids is deprecated. " + "Use orthonormal_basis_for_space instead. " + "This function will go away in 2022.", + DeprecationWarning, stacklevel=2) + + 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 = get_node_tuples(shape, n) - return mode_ids, tuple(partial(pkdo_2d, order) for order in mode_ids) - elif dims == 3: - mode_ids = get_node_tuples(shape, n) - 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): @@ -512,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 @@ -536,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 @@ -564,10 +571,14 @@ def simplex_monomial_basis_with_mode_ids(dims, n): .. versionadded:: 2018.1 """ - from modepy.shapes import get_node_tuples - mode_ids = get_node_tuples(Simplex(dims), n) + 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)) - from functools import partial return mode_ids, tuple(partial(monomial, order) for order in mode_ids) @@ -601,27 +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): - return get_basis(Simplex(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): - return get_grad_basis(Simplex(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 @@ -635,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) @@ -647,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 @@ -656,8 +687,13 @@ def tensor_product_basis(dims, basis_1d): .. versionadded:: 2017.1 """ - from modepy.shapes import Hypercube, get_node_tuples - mode_ids = get_node_tuples(Hypercube(dims), len(basis_1d)) + 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]) @@ -673,9 +709,14 @@ def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d): .. versionadded:: 2020.2 """ + 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.shapes import Hypercube, get_node_tuples - mode_ids = get_node_tuples(Hypercube(dims), len(basis_1d)) + 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( @@ -687,13 +728,21 @@ def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d): 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 @@ -736,56 +785,270 @@ def symbolicize_function(f, dims, ref_coord_var_name="r"): # }}} -# {{{ shape basis functions +# {{{ basis interface -# {{{ simplex +class BasisNotOrthonormal(Exception): + pass -@get_basis.register(Simplex) -def _(shape: Simplex, order: int): - if shape.dims <= 3: - return simplex_onb(shape.dims, order) - else: - return simplex_monomial_basis(shape.dims, order) +class Basis: + """ + .. automethod:: orthonormality_weight(r) + .. autoattribute:: mode_ids + .. autoattribute:: functions + .. autoattribute:: gradients + """ -@get_grad_basis.register(Simplex) -def _(shape: Simplex, order: int): - if shape.dims <= 3: - return grad_simplex_onb(shape.dims, order) - else: - return grad_simplex_monomial_basis(shape.dims, order) + 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__) -@get_basis_with_mode_ids.register(Simplex) -def _(shape: Simplex, order: int): - if shape.dims <= 3: - return simplex_onb_with_mode_ids(shape.dims, order) - else: - return simplex_monomial_basis_with_mode_ids(shape.dims, order) + +@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__) # }}} -# {{{ hypercube +def zerod_basis(x): + assert len(x) == 0 + x_sub = np.ones(x.shape[1:], x.dtype) + return 1 + x_sub + + +# {{{ PN bases -@get_basis.register(Hypercube) -def _(shape: Hypercube, order: int): - return legendre_tensor_product_basis(shape.dims, order) +def _pkdo_1d(order, r): + i, = order + r0, = r + return jacobi(0, 0, i, r0) -@get_grad_basis.register(Hypercube) -def _(shape: Hypercube, order: int): - return grad_legendre_tensor_product_basis(shape.dims, order) +def _grad_pkdo_1d(order, r): + i, = order + r0, = r + return (grad_jacobi(0, 0, i, r0),) -@get_basis_with_mode_ids.register(Hypercube) -def _(shape: Hypercube, order: int): - from modepy.shapes import get_node_tuples - mode_ids = get_node_tuples(shape, order) - return mode_ids, get_basis(shape, order) +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 e9cb216c..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,11 +52,16 @@ THE SOFTWARE. """ +# }}} + +from typing import List, Tuple import numpy as np import numpy.linalg as la -from modepy.shapes import Simplex, Hypercube -from modepy.shapes import get_node_count, get_node_tuples, get_unit_nodes +from functools import singledispatch, partial + +from modepy.shapes import Shape, Simplex, Hypercube +from modepy.spaces import FunctionSpace, PN, QN # {{{ equidistant nodes @@ -38,15 +74,15 @@ 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`) """ - shape = Simplex(dims) + space = PN(dims, n) if node_tuples is None: - node_tuples = get_node_tuples(shape, n) + node_tuples = node_tuples_for_space(space) else: - if len(node_tuples) != get_node_count(shape, n): + if len(node_tuples) != space.space_dim: raise ValueError("'node_tuples' list does not have the correct length") # shape: (dims, nnodes) @@ -68,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 @@ -125,11 +161,11 @@ def warp_and_blend_nodes_2d(n, node_tuples=None): except IndexError: alpha = 5/3 - shape = Simplex(2) + space = PN(2, n) if node_tuples is None: - node_tuples = get_node_tuples(shape, n) + node_tuples = node_tuples_for_space(space) else: - if len(node_tuples) != get_node_count(shape, n): + if len(node_tuples) != space.space_dim: raise ValueError("'node_tuples' list does not have the correct length") # shape: (2, nnodes) @@ -161,11 +197,11 @@ def warp_and_blend_nodes_3d(n, node_tuples=None): except IndexError: alpha = 1. - shape = Simplex(3) + space = PN(3, n) if node_tuples is None: - node_tuples = get_node_tuples(shape, n) + node_tuples = node_tuples_for_space(space) else: - if len(node_tuples) != get_node_count(shape, n): + if len(node_tuples) != space.space_dim: raise ValueError("'node_tuples' list does not have the correct length") # shape: (3, nnodes) @@ -295,20 +331,18 @@ def warp_and_blend_nodes(dims, n, node_tuples=None): 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 - """ - if dims == 0: - # NOTE: using this to maintain consistency in the 0d case - return warp_and_blend_nodes(dims, 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) @@ -320,65 +354,121 @@ def legendre_gauss_lobatto_tensor_product_nodes(dims, n): # }}} -# {{{ shape nodes +# {{{ space-based interface -# {{{ simplex +@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)) -@get_node_count.register(Simplex) -def _(shape: Simplex, order: int): - try: - from math import comb # comb is v3.8+ - node_count = comb(order + shape.dims, shape.dims) - except ImportError: - from functools import reduce - from operator import mul - node_count = reduce(mul, range(order + 1, order + shape.dims + 1), 1) \ - // reduce(mul, range(1, shape.dims + 1), 1) - return node_count +@singledispatch +def edge_clustered_nodes_for_space(space: FunctionSpace, shape: Shape): + raise NotImplementedError((type(space).__name__, type(shape).__name)) -@get_node_tuples.register(Simplex) -def _(shape: Simplex, order: int): +@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 - if shape.dims == 0: - return ((0,),) - else: - return tuple(gnitsam(order, shape.dims)) + return tuple(gnitsam(space.order, space.spatial_dim)) -@get_unit_nodes.register(Simplex) -def _(shape: Simplex, order: int): - import modepy as mp - return mp.warp_and_blend_nodes(shape.dims, order) +@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) -# {{{ hypercube +@random_nodes_for_shape.register(Simplex) +def _(shape: Simplex, nnodes: int, rng=None): + if rng is None: + rng = np.random -@get_node_count.register(Hypercube) -def _(shape: Hypercube, order: int): - return (order + 1)**shape.dims + 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 -@get_node_tuples.register(Hypercube) -def _(shape: Hypercube, order: int): +# }}} + + +# {{{ QN + +@node_tuples_for_space.register(QN) +def _(space: QN): from pytools import \ generate_nonnegative_integer_tuples_below as gnitb - if shape.dims == 0: - return ((0,),) - else: - return tuple(gnitb(order, shape.dims)) + return tuple(gnitb(space.order, space.spatial_dim)) -@get_unit_nodes.register(Hypercube) -def _(shape: Hypercube, order: int): - import modepy as mp - return mp.legendre_gauss_lobatto_tensor_product_nodes(shape.dims, order) +@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)) # }}} diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py index 4f213846..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): @@ -112,3 +132,54 @@ 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/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 index 7de126f5..5cf162b8 100644 --- a/modepy/shapes.py +++ b/modepy/shapes.py @@ -1,11 +1,15 @@ +# {{{ docstring + r""" :mod:`modepy.shapes` provides a generic description of the supported shapes (i.e. reference elements). +.. currentmodule:: modepy + .. autoclass:: Shape -.. autofunction:: get_unit_vertices -.. autofunction:: get_face_vertex_indices -.. autofunction:: get_face_map +.. autoclass:: Face +.. autofunction:: biunit_vertices_for_shape +.. autofunction:: faces_for_shape Simplices ^^^^^^^^^ @@ -17,7 +21,7 @@ Coordinates on the triangle --------------------------- -Unit coordinates :math:`(r, s)`:: +Bi-unit coordinates :math:`(r, s)` (also called 'unit' coordinates):: ^ s | @@ -29,7 +33,7 @@ | \ A-----B--> r -Vertices in unit coordinates:: +Vertices in bi-unit coordinates:: O = ( 0, 0) A = (-1, -1) @@ -58,7 +62,7 @@ Coordinates on the tetrahedron ------------------------------ -Unit coordinates :math:`(r, s, t)`:: +Bi-unit coordinates :math:`(r, s, t)` (also called 'unit' coordinates):: ^ s | @@ -75,7 +79,7 @@ (squint, and it might start making sense...) -Vertices in unit coordinates :math:`(r, s, t)`:: +Vertices in bi-unit coordinates :math:`(r, s, t)`:: O = ( 0, 0, 0) A = (-1, -1, -1) @@ -101,7 +105,7 @@ Coordinates on the square ------------------------- -Unit coordinates on :math:`(r, s)`:: +Bi-unit coordinates on :math:`(r, s)` (also called 'unit' coordinates):: ^ s | @@ -114,7 +118,7 @@ A---------B --> r -Vertices in unit coordinates:: +Vertices in bi-unit coordinates:: O = ( 0, 0) A = (-1, -1) @@ -127,12 +131,12 @@ Coordinates on the cube ----------------------- -Unit coordinates on :math:`(r, s, t)`:: +Bi-unit coordinates on :math:`(r, s, t)` (also called 'unit' coordinates):: t ^ | - B----------D + E----------G |\ |\ | \ | \ | \ | \ @@ -142,26 +146,53 @@ \ | \ | \ | \ | \| \| - E----------G + B----------D \ v r -Verties in unit coordinates:: +Vertices in bi-unit coordinates:: O = ( 0, 0, 0) A = (-1, -1, -1) - B = (-1, -1, 1) + B = ( 1, -1, -1) C = (-1, 1, -1) - D = (-1, 1, 1) - E = ( 1, -1, -1) + D = ( 1, 1, -1) + E = (-1, -1, 1) F = ( 1, -1, 1) - G = ( 1, 1, -1) + G = (-1, 1, 1) H = ( 1, 1, 1) The order of the vertices in the hypercubes follows binary counting -in ``rst``. For example, in 3D, ``A, B, C, D, ...`` is ``000, 001, 010, 011, ...``. +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 @@ -189,8 +220,9 @@ """ import numpy as np +from typing import Tuple, Callable -from functools import singledispatch +from functools import singledispatch, partial from dataclasses import dataclass @@ -199,76 +231,56 @@ @dataclass(frozen=True) class Shape: """ - .. attribute :: dims + .. attribute :: dim .. attribute :: nfaces .. attribute :: nvertices """ - dims: int - - -@singledispatch -def get_unit_vertices(shape: Shape): - """ - :returns: an :class:`~numpy.ndarray` of shape `(nvertices, dims)`. - """ - raise NotImplementedError(type(shape).__name__) - - -@singledispatch -def get_face_vertex_indices(shape: Shape): - """ - :results: a tuple of the length :attr:`Shape.nfaces`, where each entry - is a tuple of indices into the vertices returned by - :func:`get_unit_vertices`. - """ - raise NotImplementedError(type(shape).__name__) - - -@singledispatch -def get_face_map(shape: Shape, face_vertices: np.ndarray): - """ - :returns: a :class:`~collections.abc.Callable` that takes an array of - size `(dims, nnodes)` of unit nodes on the face represented by - *face_vertices* and maps them to the volume. - """ - raise NotImplementedError(type(shape).__name__) + dim: int @singledispatch -def get_quadrature(shape: Shape, order: int): +def biunit_vertices_for_shape(shape: Shape): """ - :returns: a :class:`~modepy.Quadrature` that is exact up to :math:`2 N + 1`. + :returns: an :class:`~numpy.ndarray` of shape `(dim, nvertices)`. """ raise NotImplementedError(type(shape).__name__) -@singledispatch -def get_node_count(shape: Shape, order: int): - 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 -@singledispatch -def get_node_tuples(shape: Shape, order: int): - raise NotImplementedError(type(shape).__name__) + The volume :class:`Shape` from which this face descends. + .. attribute:: face_index -@singledispatch -def get_unit_nodes(shape: Shape, order: int): - raise NotImplementedError(type(shape).__name__) + The face index in :attr:`volume_shape` of this face. + .. attribute:: volume_vertex_indices -@singledispatch -def get_basis(shape: Shape, order: int): - raise NotImplementedError(type(shape).__name__) + A tuple of indices into the vertices returned by + :func:`biunit_vertices_for_shape` for the :attr:`volume_shape`. + .. attribute:: map_to_volume -@singledispatch -def get_grad_basis(shape: Shape, order: int): - raise NotImplementedError(type(shape).__name__) + 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 get_basis_with_mode_ids(shape: Shape, order: int): +def faces_for_shape(shape: Shape): + """ + :results: a tuple of :class:`Face` representing the faces of *shape*. + """ raise NotImplementedError(type(shape).__name__) # }}} @@ -279,51 +291,57 @@ def get_basis_with_mode_ids(shape: Shape, order: int): class Simplex(Shape): @property def nfaces(self): - return self.dims + 1 + return self.dim + 1 @property def nvertices(self): return self.dim + 1 -@get_unit_vertices.register(Simplex) -def _(shape: Simplex): - from modepy.tools import unit_vertices - return unit_vertices(shape.dims) +@dataclass(frozen=True) +class _SimplexFace(Simplex, Face): + pass -@get_face_vertex_indices.register(Simplex) +@biunit_vertices_for_shape.register(Simplex) def _(shape: Simplex): - fvi = np.empty((shape.dims + 1, shape.dims), dtype=np.int) - indices = np.arange(shape.dims + 1) + result = np.empty((shape.dim, shape.dim+1), np.float64) + result.fill(-1) - for iface in range(shape.nfaces): - fvi[iface, :] = np.hstack([indices[:iface], indices[iface + 1:]]) + for i in range(shape.dim): + result[i, i+1] = 1 - return fvi + return result -@get_face_map.register(Simplex) -def _(shape: Simplex, face_vertices: np.ndarray): - dims, npoints = face_vertices.shape - if npoints != dims: +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 lambda p: origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2) + return origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2) -@get_quadrature.register(Simplex) -def _(shape: Simplex, order: int): - import modepy as mp - try: - quad = mp.XiaoGimbutasSimplexQuadrature(2*order + 1, shape.dims) - except (mp.QuadratureRuleUnavailable, ValueError): - quad = mp.GrundmannMoellerSimplexQuadrature(order, shape.dims) +@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) - return quad + 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)] # }}} @@ -333,23 +351,43 @@ def _(shape: Simplex, order: int): class Hypercube(Shape): @property def nfaces(self): - return 2 * self.dims + return 2 * self.dim @property def nvertices(self): - return 2**self.dims + return 2**self.dim + + +@dataclass(frozen=True) +class _HypercubeFace(Hypercube, Face): + pass -@get_unit_vertices.register(Hypercube) +@biunit_vertices_for_shape.register(Hypercube) def _(shape: Hypercube): from modepy.nodes import tensor_product_nodes - return tensor_product_nodes(shape.dims, np.array([-1.0, 1.0])).T + 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) -@get_face_vertex_indices.register(Hypercube) +@faces_for_shape.register(Hypercube) def _(shape: Hypercube): # FIXME: replace by nicer n-dimensional formula - return { + face_vertex_indices = { 1: ((0b0,), (0b1,)), 2: ((0b00, 0b01), (0b10, 0b11), (0b00, 0b10), (0b01, 0b11)), 3: ( @@ -362,32 +400,148 @@ def _(shape: Hypercube): (0b000, 0b001, 0b100, 0b101,), (0b010, 0b011, 0b110, 0b111,), ) - }[shape.dims] + }[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)] -@get_face_map.register(Hypercube) -def _(shape: Hypercube, face_vertices: np.ndarray): - dims, npoints = face_vertices.shape - if npoints != 2**(dims - 1): - raise ValueError("'face_vertices' has wrong shape") +# }}} - origin = face_vertices[:, 0].reshape(-1, 1) - face_basis = face_vertices[:, -2:0:-1] - origin - return lambda p: origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2) +# {{{ 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__) -@get_quadrature.register(Hypercube) -def _(shape: Hypercube, order: int): - import modepy as mp - if shape.dims == 0: - quad = mp.Quadrature(np.empty((0, 1)), np.empty((0, 1))) +@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: - from modepy.quadrature import LegendreGaussTensorProductQuadrature - quad = LegendreGaussTensorProductQuadrature(order, shape.dims) + 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) - return quad + 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 79acb876..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: @@ -170,30 +189,21 @@ def equilateral_to_unit(equi): def unit_vertices(dim): - result = np.empty((dim+1, dim), np.float64) - result.fill(-1) - - for i in range(dim): - result[i+1, i] = 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) - 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): @@ -233,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 @@ -262,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) @@ -364,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 + return shp.submesh_for_shape(shp.Hypercube(len(node_tuples[0])), node_tuples) - 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 +# }}} +# {{{ plotting helpers def plot_element_values(n, nodes, values, resample_n=None, node_tuples=None, show_nodes=False): dims = len(nodes) @@ -431,15 +333,20 @@ def plot_element_values(n, nodes, values, resample_n=None, def _evaluate_lebesgue_function(n, nodes, shape): huge_n = 30*n - from modepy.shapes import get_basis, get_node_tuples - basis = get_basis(shape, n) - equi_node_tuples = get_node_tuples(shape, huge_n) + 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 = 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) @@ -453,7 +360,7 @@ def estimate_lebesgue_constant(n, nodes, shape=None, visualize=False): `_ 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 shape: a :class:`~modepy.shapes.Shape`. :arg visualize: visualize the function that gives rise to the @@ -471,13 +378,16 @@ def estimate_lebesgue_constant(n, nodes, shape=None, visualize=False): Renamed *domain* to *shape*. """ - dims = len(nodes) + 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(dims) + shape = Simplex(dim) else: - if shape.dims != dims: - raise ValueError(f"expected {shape.dims}-dimensional nodes") + if shape.dim != dim: + raise ValueError(f"expected {shape.dim}-dimensional nodes") lebesgue_worst, equi_node_tuples, equi_nodes = \ _evaluate_lebesgue_function(n, nodes, shape) @@ -486,16 +396,9 @@ def estimate_lebesgue_constant(n, nodes, shape=None, visualize=False): if not visualize: return lebesgue_constant - if shape.dims == 2: + if shape.dim == 2: print(f"Lebesgue constant: {lebesgue_constant}") - - from modepy.shapes import Simplex, Hypercube - if isinstance(shape, Simplex): - triangles = simplex_submesh(equi_node_tuples) - elif isinstance(shape, 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 @@ -530,7 +433,7 @@ def estimate_lebesgue_constant(n, nodes, shape=None, visualize=False): ax.set_aspect("equal") plt.show() else: - raise ValueError(f"visualization is not supported in {shape.dims}D") + raise ValueError(f"visualization is not supported in {shape.dim}D") return lebesgue_constant 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_tools.py b/test/test_tools.py index fb096474..b1921171 100644 --- a/test/test_tools.py +++ b/test/test_tools.py @@ -23,7 +23,6 @@ import numpy as np import numpy.linalg as la import modepy as mp -from modepy.shapes import Simplex, Hypercube from functools import partial import pytest @@ -89,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) @@ -121,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) @@ -142,27 +142,29 @@ def estimate_resid(inner_n): # {{{ test_resampling_matrix @pytest.mark.parametrize("dims", [1, 2, 3]) -@pytest.mark.parametrize("shape_cls", [Simplex, Hypercube]) +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) def test_resampling_matrix(dims, shape_cls, ncoarse=5, nfine=10): - from modepy.shapes import get_unit_nodes, get_basis shape = shape_cls(dims) - coarse_nodes = get_unit_nodes(shape, ncoarse) - coarse_basis = get_basis(shape, ncoarse) + coarse_space = mp.space_for_shape(shape, ncoarse) + fine_space = mp.space_for_shape(shape, nfine) - fine_nodes = get_unit_nodes(shape, nfine) - fine_basis = get_basis(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 @@ -173,16 +175,14 @@ def test_resampling_matrix(dims, shape_cls, ncoarse=5, nfine=10): # {{{ test_diff_matrix @pytest.mark.parametrize("dims", [1, 2, 3]) -@pytest.mark.parametrize("shape_cls", [Simplex, Hypercube]) +@pytest.mark.parametrize("shape_cls", [mp.Simplex, mp.Hypercube]) def test_diff_matrix(dims, shape_cls, order=5): - from modepy.shapes import get_unit_nodes, get_basis, get_grad_basis 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) - nodes = get_unit_nodes(shape, order) - basis = get_basis(shape, order) - grad_basis = get_grad_basis(shape, order) - - diff_mat = mp.differentiation_matrices(basis, grad_basis, nodes) + diff_mat = mp.differentiation_matrices(basis.functions, basis.gradients, nodes) if isinstance(diff_mat, tuple): diff_mat = diff_mat[0] @@ -199,15 +199,16 @@ def test_diff_matrix(dims, shape_cls, order=5): @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) @@ -219,31 +220,29 @@ def test_diff_matrix_permutation(dims): # }}} -# {{{ test_face_mass_matrix +# {{{ face mass matrices (deprecated) @pytest.mark.parametrize("dims", [2, 3]) -@pytest.mark.parametrize("shape_cls", [Simplex, Hypercube]) -def test_modal_face_mass_matrix(dims, shape_cls, order=3): - np.set_printoptions(linewidth=200) - shape = shape_cls(dims) - - from modepy.shapes import get_unit_vertices, get_basis - vertices = get_unit_vertices(shape).T - basis = get_basis(shape, order - 1) +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) - from modepy.shapes import get_face_vertex_indices - fvi = get_face_vertex_indices(shape) + 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(shape.nfaces): - face_vertices = vertices[:, fvi[iface]] + for face in mp.faces_for_shape(shape): + face_vertices = vertices[:, face.volume_vertex_indices] - fmm = modal_face_mass_matrix(basis, order, face_vertices, shape=shape) - fmm2 = modal_face_mass_matrix(basis, order+1, face_vertices, shape=shape) + fmm = modal_face_mass_matrix( + basis.functions, order, face_vertices) + fmm2 = modal_face_mass_matrix( + basis.functions, order+1, face_vertices) 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 {iface}" + 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) @@ -252,44 +251,110 @@ def test_modal_face_mass_matrix(dims, shape_cls, order=3): @pytest.mark.parametrize("dims", [2, 3]) -@pytest.mark.parametrize("shape_cls", [Simplex, Hypercube]) -def test_nodal_face_mass_matrix(dims, shape_cls, order=3): - np.set_printoptions(linewidth=200) - volume = shape_cls(dims) - face = shape_cls(dims - 1) +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) - from modepy.shapes import get_unit_vertices, get_unit_nodes, get_basis - vertices = get_unit_vertices(volume).T - volume_nodes = get_unit_nodes(volume, order) - volume_basis = get_basis(volume, order) - face_nodes = get_unit_nodes(face, order) - - from modepy.shapes import get_face_vertex_indices - fvi = get_face_vertex_indices(volume) + 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 - for iface in range(volume.nfaces): - face_vertices = vertices[:, fvi[iface]] + 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, volume_nodes, face_nodes, order, face_vertices, - shape=volume) + volume_basis.functions, volume_nodes, + face_nodes, order, face_vertices) fmm2 = nodal_face_mass_matrix( - volume_basis, volume_nodes, face_nodes, order+1, face_vertices, - shape=volume) + 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}" + + 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( + mp.basis_for_space(face_space).functions, + mp.edge_clustered_nodes_for_space(face_space, face))) + +# }}} + + +# {{{ 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) + + +@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 < 1e-11, f"error {error:.5e} on face {iface}" + 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( - get_basis(face, order), - get_unit_nodes(face, order))) + logger.info("mass matrix:\n%s", + mp.mass_matrix(face_basis.functions, face_nodes)) # }}} @@ -298,13 +363,13 @@ def test_nodal_face_mass_matrix(dims, shape_cls, order=3): @pytest.mark.parametrize("dims", [1, 2]) @pytest.mark.parametrize("order", [3, 5, 8]) -@pytest.mark.parametrize("shape_cls", [Simplex, Hypercube]) +@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) - from modepy.shapes import get_unit_nodes - nodes = get_unit_nodes(shape, order) + nodes = mp.edge_clustered_nodes_for_space(space, shape) from modepy.tools import estimate_lebesgue_constant lebesgue_constant = estimate_lebesgue_constant(order, nodes, shape=shape) @@ -346,8 +411,8 @@ def test_estimate_lebesgue_constant(dims, order, shape_cls, 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)) @@ -356,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)