diff --git a/modepy/__init__.py b/modepy/__init__.py index c591842..c459489 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -32,12 +32,13 @@ inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face, - modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, multi_vandermonde, nodal_mass_matrix_for_face, nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, + nodal_quadrature_bilinear_form_matrix, + nodal_quadrature_test_matrix, resampling_matrix, spectral_diag_nodal_mass_matrix, vandermonde, @@ -158,13 +159,14 @@ "legendre_gauss_tensor_product_nodes", "mass_matrix", "modal_mass_matrix_for_face", - "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", "monomial_basis_for_space", "multi_vandermonde", "nodal_mass_matrix_for_face", "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", + "nodal_quadrature_bilinear_form_matrix", + "nodal_quadrature_test_matrix", "node_tuples_for_space", "orthonormal_basis_for_space", "quadrature_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 05ea349..2d3da11 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -25,6 +25,7 @@ from collections.abc import Callable, Sequence +from typing import TypeAlias from warnings import warn import numpy as np @@ -61,11 +62,9 @@ .. autofunction:: inverse_mass_matrix .. autofunction:: mass_matrix -.. autofunction:: modal_mass_matrix_for_face -.. autofunction:: nodal_mass_matrix_for_face -.. autofunction:: modal_quad_mass_matrix_for_face -.. autofunction:: nodal_quad_mass_matrix_for_face .. autofunction:: spectral_diag_nodal_mass_matrix +.. autofunction:: nodal_quadrature_test_matrix +.. autofunction:: nodal_quadrature_bilinear_form_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 @@ -78,8 +77,11 @@ """ +NodalFunction: TypeAlias = Callable[[np.ndarray], np.ndarray] + + def vandermonde( - functions: Sequence[Callable[[np.ndarray], np.ndarray]], + functions: Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: """Return a (generalized) Vandermonde matrix. @@ -163,7 +165,7 @@ def multi_vandermonde( def resampling_matrix( - basis: Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Sequence[NodalFunction], new_nodes: np.ndarray, old_nodes: np.ndarray, least_squares_ok: bool = False @@ -218,7 +220,7 @@ def is_square(m): def differentiation_matrices( - basis: Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Sequence[NodalFunction], grad_basis: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]], nodes: np.ndarray, from_nodes: np.ndarray | None = None @@ -298,7 +300,7 @@ def diff_matrix_permutation( def inverse_mass_matrix( - basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Basis | Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: """Return a matrix :math:`A=M^{-1}`, which is the inverse of the one returned @@ -317,7 +319,7 @@ def inverse_mass_matrix( "inverse mass matrix of non-orthogonal basis" ) from None - basis_functions: Sequence[Callable[[np.ndarray], np.ndarray]] = ( + basis_functions: Sequence[NodalFunction] = ( basis.functions) else: basis_functions = basis @@ -333,7 +335,7 @@ def inverse_mass_matrix( def mass_matrix( - basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Basis | Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: r"""Return a mass matrix :math:`M`, which obeys @@ -351,61 +353,106 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) -def modal_quad_mass_matrix( - quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - ) -> np.ndarray: - r"""Using the *quadrature*, provide a matrix :math:`M` that - satisfies: +def nodal_quadrature_test_matrix( + quadrature: Quadrature, + test_functions: Sequence[NodalFunction], + nodal_interp_functions: Sequence[NodalFunction], + nodes: np.ndarray | None = None, + test_function_node_map: NodalFunction | None = None + ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: .. math:: - \displaystyle (M \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, + \displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, - where :math:`\phi_i` are the *test_functions* at the nodes - :math:`r_j` of *quadrature*, with corresponding weights :math:`w_j`. + where :math:`\phi_i` are the Lagrange basis functions obtained from + *test_functions* at *nodes*, :math:`w_j` and :math:`r_j` are the weights + and nodes from *quadrature*, and :math:`u_j` are trial solution point values + at *quadrature* nodes. - .. versionadded :: 2024.2 + *test_function_node_map* is an optional argument used, for example, to map + nodes on element faces to the element volume. This does not constitute a + change in the domain of integration. This is only used to map the nodes + passed to *test_functions*. """ - modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) + if nodes is None: + nodes = quadrature.nodes - for i, test_f in enumerate(test_functions): - modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights + if len(nodal_interp_functions) != nodes.shape[1]: + raise ValueError("nodes not unisolvent with nodal_interp_functions") - return modal_mass_matrix + vdm = vandermonde(test_functions, nodes) + test_nodes = ( + test_function_node_map(quadrature.nodes) + if test_function_node_map is not None else quadrature.nodes + ) -def nodal_quad_mass_matrix( + modal_mat = vandermonde(test_functions, test_nodes).T*quadrature.weights + + return la.solve(vdm.T, modal_mat) + + +def nodal_quadrature_bilinear_form_matrix( quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - nodes: np.ndarray | None = None, + test_functions: Sequence[NodalFunction], + trial_functions: Sequence[NodalFunction], + nodal_interp_functions_test: Sequence[NodalFunction], + nodal_interp_functions_trial: Sequence[NodalFunction], + input_nodes: np.ndarray, + output_nodes: np.ndarray | None = None, + test_function_node_map: NodalFunction | None = None, ) -> np.ndarray: - r"""Using the *quadrature*, provide a matrix :math:`M` that - satisfies: + r"""Using *quadrature*, provide a matrix :math:`A` defined as: .. math:: - \displaystyle (M \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, + \displaystyle A_{ij} = \sum_k \psi_i(r_k) \phi_j(r_k) w_k, - where :math:`\phi_i` are the (volume) Lagrange basis functions obtained from - *test_functions* at *nodes*, :math:`w_i` and :math:`r_i` are the - weights and nodes from *qudarature*, and :math:`u_j` are the point - values of the trial function at the nodes of *quadrature*. + where :math:`\psi_i` are Lagrange basis functions obtained from + *test_functions*, :math:`\phi_j` are Lagrange basis functions obtained from + *trial_functions*, :math:`r_k` and :math:`w_k` are nodes and weights from + *quadrature*. The matrix :math:`A` satisfies - If *nodes* is not provided, use *quadrature*'s nodes. + .. math:: - .. versionadded :: 2024.2 + \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, + + where :math:`u_i` are nodal coefficients. + + *test_function_node_map* is an optional argument used, for example, to map + nodes on element faces to the element volume. This does not constitute a + change in the domain of integration. This is only used to map the nodes + passed to *test_functions*. """ - if nodes is None: - nodes = quadrature.nodes + if output_nodes is None: + output_nodes = input_nodes - if len(test_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with test_functions") + if len(nodal_interp_functions_test) != output_nodes.shape[1]: + raise ValueError( + "output_nodes not unisolvent with nodal_test_interp_functions") - vdm = vandermonde(test_functions, nodes) + if len(nodal_interp_functions_trial) != input_nodes.shape[1]: + raise ValueError( + "input_nodes not unisolvent with nodal_trial_interp_functions") - return la.solve(vdm.T, - modal_quad_mass_matrix(quadrature, test_functions)) + mapped_nodes = ( + test_function_node_map(quadrature.nodes) + if test_function_node_map is not None else quadrature.nodes + ) + + modal_operator = np.einsum( + "qi,qj,q->ij", + vandermonde(test_functions, mapped_nodes), + vandermonde(trial_functions, quadrature.nodes), + quadrature.weights + ) + + input_vdm = vandermonde(nodal_interp_functions_trial, input_nodes) + output_vdm = vandermonde(nodal_interp_functions_test, output_nodes) + + return la.solve(output_vdm.T, modal_operator @ la.inv(input_vdm)) def spectral_diag_nodal_mass_matrix( @@ -429,25 +476,51 @@ def spectral_diag_nodal_mass_matrix( return quadrature.weights -def modal_mass_matrix_for_face( - face: Face, face_quad: Quadrature, - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]] +# {{{ deprecated remove in 2026-ish + +def modal_quad_mass_matrix( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: - r"""Using the quadrature *face_quad*, provide a matrix :math:`M_f` that - satisfies: + from warnings import warn + warn("`modal_quad_mass_matrix` is deprecated and will stop working in " + "2026.", stacklevel=1) + modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) - .. math:: + for i, test_f in enumerate(test_functions): + modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights - \displaystyle {(M_f)}_{ij} \approx \int_F \phi_i(r) \psi_j(r) dr, + return modal_mass_matrix - where :math:`\phi_i` are the (volume) basis functions - *test_functions*, and :math:`\psi_j` are the (surface) - *trial_functions*. - .. versionadded:: 2020.3 - """ +def nodal_quad_mass_matrix( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + nodes: np.ndarray | None = None, + ) -> np.ndarray: + from warnings import warn + warn("`nodal_quad_mass_matrix` is deprecated and will stop working in " + "2026. Consider switching to `nodal_quad_bilinear_form`", stacklevel=1) + if nodes is None: + nodes = quadrature.nodes + if len(test_functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with test_functions") + + vdm = vandermonde(test_functions, nodes) + + return la.solve(vdm.T, + modal_quad_mass_matrix(quadrature, test_functions)) + + +def modal_mass_matrix_for_face( + face: Face, face_quad: Quadrature, + trial_functions: Sequence[NodalFunction], + test_functions: Sequence[NodalFunction] + ) -> np.ndarray: + from warnings import warn + warn("`modal_mass_matrix_for_face` is deprecated and will stop working in " + "2026.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) result = np.empty((len(test_functions), len(trial_functions))) @@ -462,23 +535,14 @@ def modal_mass_matrix_for_face( def nodal_mass_matrix_for_face( face: Face, face_quad: Quadrature, - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - volume_nodes: np.ndarray, face_nodes: np.ndarray + trial_functions: Sequence[NodalFunction], + test_functions: Sequence[NodalFunction], + volume_nodes: np.ndarray, + face_nodes: np.ndarray ) -> np.ndarray: - r"""Using the quadrature *face_quad*, provide a matrix :math:`M_f` that - satisfies: - - .. math:: - - \displaystyle {(M_f)}_{ij} \approx \int_F \phi_i(r) \psi_j(r) dr, - - where :math:`\phi_i` are the (volume) Lagrange basis functions obtained from - *test_functions* at *volume_nodes*, and :math:`\psi_j` are the (surface) - Lagrange basis functions obtained from *trial_functions* at *face_nodes*. - - .. versionadded :: 2020.3 - """ + from warnings import warn + warn("`nodal_mass_matrix_for_face` is deprecated and will stop working in " + "2026. Please use `nodal_quad_bilinear_form` instead.", stacklevel=1) face_vdm = vandermonde(trial_functions, face_nodes) vol_vdm = vandermonde(test_functions, volume_nodes) @@ -497,21 +561,11 @@ def nodal_mass_matrix_for_face( def modal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + test_functions: Sequence[NodalFunction], ) -> np.ndarray: - r"""Using the quadrature *face_quad*, provide a matrix :math:`M_f` that - satisfies: - - .. math:: - - \displaystyle (M_f \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, - - where :math:`\phi_i` are the *test_functions*, :math:`w_j` and :math:`r_j` are - the weights and nodes from *face_quad*, and :math:`u_j` are the point values of - the trial function at the nodes of *face_quad*. - - .. versionadded :: 2024.2 - """ + from warnings import warn + warn("`modal_quad_mass_matrix_for_face` is deprecated and will stop working " + "in 2025.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) vol_modal_mass_matrix = np.empty((len(test_functions), len(face_quad.weights))) @@ -524,30 +578,22 @@ def modal_quad_mass_matrix_for_face( def nodal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + test_functions: Sequence[NodalFunction], volume_nodes: np.ndarray, ) -> np.ndarray: - r"""Using the quadrature *face_quad*, provide a matrix :math:`M_f` that - satisfies: - - .. math:: - - \displaystyle (M_f \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, - - where :math:`\phi_i` are the (volume) Lagrange basis functions obtained from - *test_functions* at *volume_nodes*, :math:`w_j` and :math:`r_j` are the - weights and nodes from *face_quad*, and :math:`u_j` are the point - values of the trial function at the nodes of *face_quad*. - - .. versionadded :: 2021.2 - """ + from warnings import warn + warn("`nodal_quad_mass_matrix_for_face` is deprecated and will stop working " + "in 2025. Consider using `nodal_quad_bilinear_form` instead", + stacklevel=1) if len(test_functions) != volume_nodes.shape[1]: raise ValueError("volume_nodes not unisolvent with test_functions") vol_vdm = vandermonde(test_functions, volume_nodes) return la.solve(vol_vdm.T, - modal_quad_mass_matrix_for_face(face, face_quad, test_functions)) + modal_quad_mass_matrix_for_face(face, face_quad, + test_functions)) +# }}} # vim: foldmethod=marker diff --git a/modepy/modes.py b/modepy/modes.py index 05a3134..917da1c 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -31,6 +31,7 @@ from functools import partial, singledispatch from typing import ( TYPE_CHECKING, + TypeAlias, TypeVar, cast, ) @@ -46,6 +47,7 @@ RealValueT = TypeVar("RealValueT", np.ndarray, "pymbolic.primitives.Expression", float) + """:class:`~typing.TypeVar` for basis function inputs and outputs.""" __doc__ = """This functionality provides sets of basis functions for the @@ -745,8 +747,8 @@ def symbolicize_function( # {{{ basis interface -BasisFunctionType = Callable[[np.ndarray], np.ndarray] -BasisGradientType = Callable[[np.ndarray], tuple[np.ndarray, ...]] +BasisFunctionType: TypeAlias = Callable[[np.ndarray], np.ndarray] +BasisGradientType: TypeAlias = Callable[[np.ndarray], tuple[np.ndarray, ...]] class BasisNotOrthonormal(Exception): diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 68bd492..6e326fe 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -118,6 +118,69 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: assert np.max(err) < 1e-14 +@pytest.mark.parametrize("shape_cls", [mp.Hypercube, mp.Simplex]) +@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("order", [0, 1, 2, 4]) +@pytest.mark.parametrize("nodes_on_bdry", [False, True]) +@pytest.mark.parametrize("test_weak_d_dr", [False, True]) +def test_bilinear_forms( + shape_cls: type[mp.Shape], + dim: int, + order: int, + nodes_on_bdry: bool, + test_weak_d_dr: bool + ) -> None: + shape = shape_cls(dim) + space = mp.space_for_shape(shape, order) + + quad_space = mp.space_for_shape(shape, 2*order) + quad = mp.quadrature_for_space(quad_space, shape) + + if nodes_on_bdry: + nodes = mp.edge_clustered_nodes_for_space(space, shape) + else: + if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): + nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order) + elif isinstance(shape, mp.Simplex): + nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes + else: + raise AssertionError() + + basis = mp.orthonormal_basis_for_space(space, shape) + + if test_weak_d_dr and order not in [0, 1]: + mass_inv = mp.inverse_mass_matrix(basis, nodes) + + for ax in range(dim): + f = 1 - nodes[ax]**2 + fp = -2*nodes[ax] + + weak_operator = mp.nodal_quadrature_bilinear_form_matrix( + quadrature=quad, + test_functions=basis.derivatives(ax), + trial_functions=basis.functions, + nodal_interp_functions_test=basis.functions, + nodal_interp_functions_trial=basis.functions, + input_nodes=nodes, + ) + + err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) + assert err <= 1e-12 + else: + quad_mass_mat = mp.nodal_quadrature_bilinear_form_matrix( + quadrature=quad, + test_functions=basis.functions, + trial_functions=basis.functions, + nodal_interp_functions_test=basis.functions, + nodal_interp_functions_trial=basis.functions, + input_nodes=nodes + ) + + vdm_mass_mat = mp.mass_matrix(basis, nodes) + err = la.norm(quad_mass_mat - vdm_mass_mat) / la.norm(vdm_mass_mat) + assert err < 1e-14 + + # You can test individual routines by typing # $ python test_modes.py 'test_routine()' diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py index ef57566..2d22f8a 100644 --- a/modepy/test/test_tools.py +++ b/modepy/test/test_tools.py @@ -295,36 +295,7 @@ def test_diff_matrix_permutation(dims): @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, vol_shape) - - 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) - 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): +def test_nodal_quadrature_bilinear_form_matrix_for_face(dims, shape_cls, order=3): vol_shape = shape_cls(dims) vol_space = mp.space_for_shape(vol_shape, order) @@ -332,8 +303,8 @@ def test_nodal_mass_matrix_for_face(dims, shape_cls, order=3): volume_basis = mp.basis_for_space(vol_space, vol_shape) from modepy.matrices import ( - nodal_mass_matrix_for_face, - nodal_quad_mass_matrix_for_face, + nodal_quadrature_bilinear_form_matrix, + nodal_quadrature_test_matrix, ) for face in mp.faces_for_shape(vol_shape): face_space = mp.space_for_shape(face, order) @@ -342,11 +313,25 @@ def test_nodal_mass_matrix_for_face(dims, shape_cls, order=3): 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_quad_mass_matrix_for_face( - face, face_quad2, volume_basis.functions, volume_nodes) + + fmm = nodal_quadrature_bilinear_form_matrix( + quadrature=face_quad, + test_functions=volume_basis.functions, + trial_functions=face_basis.functions, + nodal_interp_functions_test=volume_basis.functions, + nodal_interp_functions_trial=face_basis.functions, + input_nodes=face_nodes, + output_nodes=volume_nodes, + test_function_node_map=face.map_to_volume + ) + + fmm2 = nodal_quadrature_test_matrix( + quadrature=face_quad2, + test_functions=volume_basis.functions, + nodal_interp_functions=volume_basis.functions, + nodes=volume_nodes, + test_function_node_map=face.map_to_volume + ) for f_face in face_basis.functions: fval_nodal = f_face(face_nodes)