diff --git a/doc/tools.rst b/doc/tools.rst index 1d64d3fa088faa08059aae1610a48270ba0470c1..b5b00069ca29006cfd399f0c62cb3bbb4c6eda2c 100644 --- a/doc/tools.rst +++ b/doc/tools.rst @@ -8,7 +8,7 @@ Version query .. data:: VERSION - A tuple like *(2013,1,1)* indicating modepy's version. + A tuple like *(2013, 1, 1)* indicating modepy's version. .. data:: VERSION_TEXT @@ -17,49 +17,13 @@ Version query Interpolation matrices ---------------------- -.. currentmodule:: modepy - -.. autofunction:: vandermonde - -Vandermonde matrices are very useful to concisely express interpolation. For -instance, one may use the inverse :math:`V^{-1}` of a Vandermonde matrix -:math:`V` to map nodal values (i.e. point values at the nodes corresponding to -:math:`V`) into modal (i.e. basis) coefficients. :math:`V` itself maps modal -coefficients to nodal values. This allows easy resampling: - -.. autofunction:: resampling_matrix - -Vandermonde matrices also obey the following relationship useful for obtaining -point interpolants: - -.. math:: - - V^T [\text{interpolation coefficents to point $x$}] = \phi_i(x), - -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 - -Differentiation is also convenient to express by using :math:`V^{-1}` to -obtain modal values and then using a Vandermonde matrix for the derivatives -of the basis to return to nodal values. - -.. autofunction:: differentiation_matrices +.. automodule:: modepy.matrices Modal decay/residual -------------------- .. automodule:: modepy.modal_decay -.. autofunction:: fit_modal_decay -.. autofunction:: estimate_relative_expansion_residual - Interpolation quality --------------------- diff --git a/modepy/matrices.py b/modepy/matrices.py index 4e53001a530d853ef9a31d1577089be7bfc217c8..34aaaa31314b449e63dddcfb33eeb98881b3e3a7 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -27,6 +27,44 @@ import numpy as np import numpy.linalg as la +__doc__ = r""" +.. currentmodule:: modepy + +.. autofunction:: vandermonde + +Vandermonde matrices are very useful to concisely express interpolation. For +instance, one may use the inverse :math:`V^{-1}` of a Vandermonde matrix +:math:`V` to map nodal values (i.e. point values at the nodes corresponding to +:math:`V`) into modal (i.e. basis) coefficients. :math:`V` itself maps modal +coefficients to nodal values. This allows easy resampling: + +.. autofunction:: resampling_matrix + +Vandermonde matrices also obey the following relationship useful for obtaining +point interpolants: + +.. math:: + + V^T [\text{interpolation coefficents to point $x$}] = \phi_i(x), + +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 + +Differentiation is also convenient to express by using :math:`V^{-1}` to +obtain modal values and then using a Vandermonde matrix for the derivatives +of the basis to return to nodal values. + +.. autofunction:: differentiation_matrices +""" + + def vandermonde(functions, nodes): """Return a (generalized) Vandermonde matrix. diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py index 9004123baf80ce42e25ce079826193c43e1cae5e..dee45fee15614f4c461be570f4d9c38b4b3eeb11 100644 --- a/modepy/modal_decay.py +++ b/modepy/modal_decay.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2010-2012 Andreas Kloeckner" @@ -24,6 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range + import numpy as np import numpy.linalg as la @@ -39,9 +39,42 @@ The method implemented in this module follows this article: http://arxiv.org/abs/1102.3190 .. versionadded:: 2013.2 + +.. autofunction:: simplex_interp_error_coefficient_estimator_matrix +.. autofunction:: fit_modal_decay +.. autofunction:: estimate_relative_expansion_residual """ +def simplex_interp_error_coefficient_estimator_matrix( + unit_nodes, order, n_tail_orders): + """Return a matrix :math:`C` that, when multiplied by a vector of nodal values, + yields the coeffiicients belonging to the basis functions of the *n_tail_orders* + highest orders. + + The 2-norm of the resulting coefficents can be used as an estimate of the + interpolation error. + + .. versionadded:: 2018.1 + """ + + from modepy.matrices import vandermonde + from modepy.modes import simplex_onb_with_mode_ids + + dim, nunit_nodes = unit_nodes.shape + + mode_ids, basis = simplex_onb_with_mode_ids(dim, order) + vdm = vandermonde(basis, unit_nodes) + vdm_inv = la.inv(vdm) + + order_vector = np.array([sum(mode_id) for mode_id in mode_ids]) + + max_order = np.max(order_vector) + assert max_order == order + + return vdm_inv[order_vector > max_order-n_tail_orders] + + def make_mode_number_vector(mode_order_tuples, ignored_modes): node_cnt = len(mode_order_tuples) diff --git a/modepy/modes.py b/modepy/modes.py index e6dfcc197777e153a5dfe2d1c658db9f21a8f864..f8416298a2691cc8242f5ac7ebb9109e59325de3 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -61,10 +61,14 @@ Dimension-independent basis getters for simplices 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 @@ -438,15 +442,17 @@ def grad_monomial(order, rst): # {{{ dimension-independent interface for simplices -def simplex_onb(dims, n): +def simplex_onb_with_mode_ids(dims, n): """Return a list of orthonormal basis functions in dimension *dims* of maximal total degree *n*. - :returns: a class:`tuple` of functions, each of which - accepts arrays of shape *(dims, npts)* - and return the function values as an array of size *npts*. - 'Scalar' evaluation, by passing just one vector of length *dims*, - is also supported. + :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple` + of functions, each of which accepts arrays of shape *(dims, npts)* and + return the function values as an array of size *npts*. 'Scalar' + evaluation, by passing just one vector of length *dims*, is also supported. + *mode_ids* is a tuple of the same length as *basis*, where each entry + is a tuple of integers describing the order of the mode along the coordinate + axes. See the following publications: @@ -454,9 +460,7 @@ def simplex_onb(dims, n): * |koornwinder-ref| * |dubiner-ref| - .. versionchanged:: 2013.2 - - Made return value a tuple, to make bases hashable. + ... versionadded: 2018.1 """ from functools import partial from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ @@ -469,18 +473,45 @@ def simplex_onb(dims, n): else: return np.ones(x.shape[1]) - return (zerod_basis,) + return ((0,),), (zerod_basis,) elif dims == 1: - return tuple(partial(jacobi, 0, 0, i) for i in range(n+1)) + mode_ids = tuple(range(n+1)) + return mode_ids, tuple(partial(jacobi, 0, 0, i) for i in mode_ids) elif dims == 2: - return tuple(partial(pkdo_2d, order) for order in gnitstam(n, dims)) + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(pkdo_2d, order) for order in mode_ids) elif dims == 3: - return tuple(partial(pkdo_3d, order) for order in gnitstam(n, dims)) + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(pkdo_3d, order) for order in mode_ids) else: raise NotImplementedError("%d-dimensional bases" % dims) +def simplex_onb(dims, n): + """Return a list of orthonormal basis functions in dimension *dims* of maximal + total degree *n*. + + :returns: a class:`tuple` of functions, each of which + accepts arrays of shape *(dims, npts)* + and return the function values as an array of size *npts*. + 'Scalar' evaluation, by passing just one vector of length *dims*, + is also supported. + + See the following publications: + + * |proriol-ref| + * |koornwinder-ref| + * |dubiner-ref| + + .. versionchanged:: 2013.2 + + Made return value a tuple, to make bases hashable. + """ + mode_ids, basis = simplex_onb_with_mode_ids(dims, n) + return basis + + def grad_simplex_onb(dims, n): """Return the gradients of the functions returned by :func:`simplex_onb`. @@ -515,6 +546,27 @@ def grad_simplex_onb(dims, n): raise NotImplementedError("%d-dimensional bases" % dims) +def simplex_monomial_basis_with_mode_ids(dims, n): + """Return a list of monomial basis functions in dimension *dims* of maximal + total degree *n*. + + :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple` + of functions, each of which accepts arrays of shape *(dims, npts)* and + return the function values as an array of size *npts*. 'Scalar' + evaluation, by passing just one vector of length *dims*, is also supported. + *mode_ids* is a tuple of the same length as *basis*, where each entry + is a tuple of integers describing the order of the mode along the coordinate + axes. + + .. versionadded:: 2018.1 + """ + from functools import partial + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ + as gnitstam + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(monomial, order) for order in mode_ids) + + def simplex_monomial_basis(dims, n): """Return a list of monomial basis functions in dimension *dims* of maximal total degree *n*. @@ -527,11 +579,8 @@ def simplex_monomial_basis(dims, n): .. versionadded:: 2016.1 """ - - from functools import partial - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - return tuple(partial(monomial, order) for order in gnitstam(n, dims)) + mode_ids, basis = simplex_monomial_basis_with_mode_ids(dims, n) + return basis def grad_simplex_monomial_basis(dims, n):