From 88ac3b56eaa207086d6f0e5f252ec5bb490dd579 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 14:55:19 -0600 Subject: [PATCH 01/28] add nodal and modal bilinear form routines --- modepy/matrices.py | 66 ++++++++++++++++++++++++++++++++++ modepy/test/test_matrices.py | 68 +++++++++++++++++++++++++++++++++++- 2 files changed, 133 insertions(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 05ea3494..1bf7ce35 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -550,4 +550,70 @@ def nodal_quad_mass_matrix_for_face( modal_quad_mass_matrix_for_face(face, face_quad, test_functions)) +def modal_quad_bilinear_form( + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + quadrature: Quadrature + ) -> np.ndarray: + r"""Using the *quadrature*, provide a matrix :math:`A` that + satisfies: + + .. math:: + + \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`. + + Seeks to generalize functionality found in, e.g. + :func:`modal_quad_mass_matrix`. + """ + modal_operator = np.empty((len(test_functions), len(quadrature.weights))) + + for i, test_f in enumerate(test_functions): + modal_operator[i] = test_f(quadrature.nodes) * quadrature.weights + + return modal_operator + + +def nodal_quad_bilinear_form( + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + quadrature: Quadrature, + nodes: np.ndarray, + uses_quadrature_domain: bool = False + ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + + .. math:: + + \displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_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 point values of the trial + function at either *nodes* or the *quadrature* nodes depending on whether + a quadrature domain is used (as signified by *uses_quadrature_domain*). + + If *uses_quadrature_domain* is set to False, then an interpolation operator + is used to make the operator :math:`N \times N` where :math:`N` is the + number of nodes in *nodes*. Otherwise, the operator has shape :math:`N + \times N_q` where :math:`N_q` is the number of nodes associated with + *quadrature*. + + Seeks to generalize functionality found in, e.g. + :func:`nodal_quad_mass_matrix`. + """ + if len(test_functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with test_functions") + + vdm_out = vandermonde(trial_functions, nodes) + + modal_operator = modal_quad_bilinear_form(test_functions, quadrature) + if uses_quadrature_domain: + return la.solve(vdm_out.T, modal_operator) + else: + vdm_in = vandermonde(trial_functions, quadrature.nodes) + nodal_interp_mat = la.solve(vdm_out.T, vdm_in.T) + return la.solve(vdm_out.T, modal_operator) @ nodal_interp_mat.T + # vim: foldmethod=marker diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 68bd492f..e2585a6e 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -1,5 +1,7 @@ from __future__ import annotations +from modepy.matrices import nodal_quad_bilinear_form, nodal_quad_mass_matrix + __copyright__ = "Copyright (C) 2024 University of Illinois Board of Trustees" @@ -24,7 +26,7 @@ """ -from typing import cast +from typing import Type, cast import numpy as np import numpy.linalg as la @@ -118,6 +120,70 @@ 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_derivatives", [False, True]) +def test_bilinear_forms( + shape_cls: Type[mp.Shape], + dim: int, + order: int, + nodes_on_bdry: bool, + test_derivatives: 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 isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): + if nodes_on_bdry: + nodes = mp.legendre_gauss_lobatto_tensor_product_nodes( + shape.dim, order, + ) + else: + nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order) + elif isinstance(shape, mp.Simplex): + if nodes_on_bdry: + nodes = mp.warp_and_blend_nodes(shape.dim, order) + else: + nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes + else: + raise AssertionError() + + basis = mp.orthonormal_basis_for_space(space, shape) + + if test_derivatives 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 = nodal_quad_bilinear_form( + basis.derivatives(ax), + basis.functions, + quad, + nodes, + uses_quadrature_domain=False) + + err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) + assert err <= 1e-12 + else: + quad_mass_mat = nodal_quad_bilinear_form( + basis.functions, + basis.functions, + quad, + nodes, + uses_quadrature_domain=False) + + 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()' From ee62ab070d7844067bc222cd4ad7736deac11e54 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 15:04:42 -0600 Subject: [PATCH 02/28] use resampling routine; make ruff and mypy happy --- modepy/__init__.py | 4 ++++ modepy/matrices.py | 7 ++++--- modepy/test/test_matrices.py | 10 ++++------ 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index c591842a..8c78884f 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -32,10 +32,12 @@ inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face, + modal_quad_bilinear_form, modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, multi_vandermonde, nodal_mass_matrix_for_face, + nodal_quad_bilinear_form, nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, resampling_matrix, @@ -158,11 +160,13 @@ "legendre_gauss_tensor_product_nodes", "mass_matrix", "modal_mass_matrix_for_face", + "modal_quad_bilinear_form", "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", "monomial_basis_for_space", "multi_vandermonde", "nodal_mass_matrix_for_face", + "nodal_quad_bilinear_form", "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", "node_tuples_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 1bf7ce35..1152877f 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -66,6 +66,8 @@ .. autofunction:: modal_quad_mass_matrix_for_face .. autofunction:: nodal_quad_mass_matrix_for_face .. autofunction:: spectral_diag_nodal_mass_matrix +.. autofunction:: modal_quad_bilinear_form +.. autofunction:: nodal_quad_bilinear_form Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives @@ -612,8 +614,7 @@ def nodal_quad_bilinear_form( if uses_quadrature_domain: return la.solve(vdm_out.T, modal_operator) else: - vdm_in = vandermonde(trial_functions, quadrature.nodes) - nodal_interp_mat = la.solve(vdm_out.T, vdm_in.T) - return la.solve(vdm_out.T, modal_operator) @ nodal_interp_mat.T + resampler = resampling_matrix(trial_functions, quadrature.nodes, nodes) + return la.solve(vdm_out.T, modal_operator) @ resampler # vim: foldmethod=marker diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index e2585a6e..5413c951 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -1,7 +1,5 @@ from __future__ import annotations -from modepy.matrices import nodal_quad_bilinear_form, nodal_quad_mass_matrix - __copyright__ = "Copyright (C) 2024 University of Illinois Board of Trustees" @@ -26,7 +24,7 @@ """ -from typing import Type, cast +from typing import cast import numpy as np import numpy.linalg as la @@ -126,7 +124,7 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: @pytest.mark.parametrize("nodes_on_bdry", [False, True]) @pytest.mark.parametrize("test_derivatives", [False, True]) def test_bilinear_forms( - shape_cls: Type[mp.Shape], + shape_cls: type[mp.Shape], dim: int, order: int, nodes_on_bdry: bool, @@ -162,7 +160,7 @@ def test_bilinear_forms( f = 1 - nodes[ax]**2 fp = -2*nodes[ax] - weak_operator = nodal_quad_bilinear_form( + weak_operator = mp.nodal_quad_bilinear_form( basis.derivatives(ax), basis.functions, quad, @@ -172,7 +170,7 @@ def test_bilinear_forms( err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) assert err <= 1e-12 else: - quad_mass_mat = nodal_quad_bilinear_form( + quad_mass_mat = mp.nodal_quad_bilinear_form( basis.functions, basis.functions, quad, From 26e1258404f429af1abacdbfa3ee5ee39ebb799f Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 15:15:11 -0600 Subject: [PATCH 03/28] clean it up --- modepy/matrices.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 1152877f..4b8c0c43 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -606,15 +606,17 @@ def nodal_quad_bilinear_form( :func:`nodal_quad_mass_matrix`. """ if len(test_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with test_functions") + raise ValueError("volume_nodes not unisolvent with trial_functions") vdm_out = vandermonde(trial_functions, nodes) - modal_operator = modal_quad_bilinear_form(test_functions, quadrature) + nodal_operator = la.solve( + vdm_out.T, modal_quad_bilinear_form(test_functions, quadrature)) + if uses_quadrature_domain: - return la.solve(vdm_out.T, modal_operator) + return nodal_operator else: - resampler = resampling_matrix(trial_functions, quadrature.nodes, nodes) - return la.solve(vdm_out.T, modal_operator) @ resampler + return nodal_operator @ resampling_matrix( + trial_functions, quadrature.nodes, nodes) # vim: foldmethod=marker From c99dc720135a88ed7f5bf1011d4de67bd52b4ae1 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 15:19:20 -0600 Subject: [PATCH 04/28] test_derivatives -> test_weak_d_dr for clarity --- modepy/test/test_matrices.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 5413c951..e69d5d23 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -122,13 +122,13 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: @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_derivatives", [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_derivatives: bool + test_weak_d_dr: bool ) -> None: shape = shape_cls(dim) space = mp.space_for_shape(shape, order) @@ -153,7 +153,7 @@ def test_bilinear_forms( basis = mp.orthonormal_basis_for_space(space, shape) - if test_derivatives and order not in [0, 1]: + if test_weak_d_dr and order not in [0, 1]: mass_inv = mp.inverse_mass_matrix(basis, nodes) for ax in range(dim): From 19b1cda73233c8c927e85f8f56c5706c4d2a6fea Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 15:44:26 -0600 Subject: [PATCH 05/28] update moda/nodal quadrature mass routines to use new routines --- modepy/matrices.py | 156 ++++++++++++++++++++++----------------------- 1 file changed, 75 insertions(+), 81 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 4b8c0c43..9622bd5f 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -353,6 +353,74 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) +def modal_quad_bilinear_form( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]] + ) -> np.ndarray: + r"""Using the *quadrature*, provide a matrix :math:`A` that + satisfies: + + .. math:: + + \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`. + + Seeks to generalize functionality found in, e.g. + :func:`modal_quad_mass_matrix`. + """ + modal_operator = np.empty((len(test_functions), len(quadrature.weights))) + + for i, test_f in enumerate(test_functions): + modal_operator[i] = test_f(quadrature.nodes) * quadrature.weights + + return modal_operator + + +def nodal_quad_bilinear_form( + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + quadrature: Quadrature, + nodes: np.ndarray, + uses_quadrature_domain: bool = False + ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + + .. math:: + + \displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_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 point values of the trial + function at either *nodes* or the *quadrature* nodes depending on whether + a quadrature domain is used (as signified by *uses_quadrature_domain*). + + If *uses_quadrature_domain* is set to False, then an interpolation operator + is used to make the operator :math:`N \times N` where :math:`N` is the + number of nodes in *nodes*. Otherwise, the operator has shape :math:`N + \times N_q` where :math:`N_q` is the number of nodes associated with + *quadrature*. Default behavior is to assume a quadrature domain is not used. + + Seeks to generalize functionality found in, e.g. + :func:`nodal_quad_mass_matrix`. + """ + if len(test_functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with trial_functions") + + vdm_out = vandermonde(trial_functions, nodes) + + nodal_operator = la.solve( + vdm_out.T, modal_quad_bilinear_form(quadrature, test_functions)) + + if uses_quadrature_domain: + return nodal_operator + else: + return nodal_operator @ resampling_matrix( + trial_functions, quadrature.nodes, nodes) + + def modal_quad_mass_matrix( quadrature: Quadrature, test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], @@ -369,12 +437,7 @@ def modal_quad_mass_matrix( .. versionadded :: 2024.2 """ - modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) - - for i, test_f in enumerate(test_functions): - modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights - - return modal_mass_matrix + return modal_quad_bilinear_form(quadrature, test_functions) def nodal_quad_mass_matrix( @@ -401,13 +464,12 @@ def nodal_quad_mass_matrix( 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)) + return nodal_quad_bilinear_form( + test_functions, + test_functions, + quadrature, + nodes, + uses_quadrature_domain=True) def spectral_diag_nodal_mass_matrix( @@ -551,72 +613,4 @@ def nodal_quad_mass_matrix_for_face( return la.solve(vol_vdm.T, modal_quad_mass_matrix_for_face(face, face_quad, test_functions)) - -def modal_quad_bilinear_form( - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - quadrature: Quadrature - ) -> np.ndarray: - r"""Using the *quadrature*, provide a matrix :math:`A` that - satisfies: - - .. math:: - - \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`. - - Seeks to generalize functionality found in, e.g. - :func:`modal_quad_mass_matrix`. - """ - modal_operator = np.empty((len(test_functions), len(quadrature.weights))) - - for i, test_f in enumerate(test_functions): - modal_operator[i] = test_f(quadrature.nodes) * quadrature.weights - - return modal_operator - - -def nodal_quad_bilinear_form( - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - quadrature: Quadrature, - nodes: np.ndarray, - uses_quadrature_domain: bool = False - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: - - .. math:: - - \displaystyle (A \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_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 point values of the trial - function at either *nodes* or the *quadrature* nodes depending on whether - a quadrature domain is used (as signified by *uses_quadrature_domain*). - - If *uses_quadrature_domain* is set to False, then an interpolation operator - is used to make the operator :math:`N \times N` where :math:`N` is the - number of nodes in *nodes*. Otherwise, the operator has shape :math:`N - \times N_q` where :math:`N_q` is the number of nodes associated with - *quadrature*. - - Seeks to generalize functionality found in, e.g. - :func:`nodal_quad_mass_matrix`. - """ - if len(test_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with trial_functions") - - vdm_out = vandermonde(trial_functions, nodes) - - nodal_operator = la.solve( - vdm_out.T, modal_quad_bilinear_form(test_functions, quadrature)) - - if uses_quadrature_domain: - return nodal_operator - else: - return nodal_operator @ resampling_matrix( - trial_functions, quadrature.nodes, nodes) - # vim: foldmethod=marker From 6cbcfd5114047aea46544c75d89abc05fa1a1514 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 5 Dec 2024 15:53:35 -0600 Subject: [PATCH 06/28] fix failing docs --- modepy/matrices.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 9622bd5f..8c672f50 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -366,9 +366,6 @@ def modal_quad_bilinear_form( where :math:`\phi_i` are the *test_functions* at the nodes :math:`r_j` of *quadrature*, with corresponding weights :math:`w_j`. - - Seeks to generalize functionality found in, e.g. - :func:`modal_quad_mass_matrix`. """ modal_operator = np.empty((len(test_functions), len(quadrature.weights))) @@ -402,9 +399,6 @@ def nodal_quad_bilinear_form( number of nodes in *nodes*. Otherwise, the operator has shape :math:`N \times N_q` where :math:`N_q` is the number of nodes associated with *quadrature*. Default behavior is to assume a quadrature domain is not used. - - Seeks to generalize functionality found in, e.g. - :func:`nodal_quad_mass_matrix`. """ if len(test_functions) != nodes.shape[1]: raise ValueError("volume_nodes not unisolvent with trial_functions") From 8ea7aedaef64483964bebbffdd077ebf2ba4f6d9 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 6 Dec 2024 15:11:32 -0600 Subject: [PATCH 07/28] fix failing ruff test --- contrib/weights-from-festa-sommariva.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/contrib/weights-from-festa-sommariva.py b/contrib/weights-from-festa-sommariva.py index baab0f9e..ca720466 100644 --- a/contrib/weights-from-festa-sommariva.py +++ b/contrib/weights-from-festa-sommariva.py @@ -76,17 +76,17 @@ def generate_festa_sommariva_quadrature_rules(outfile): # # DEGREE: - _re_degree = re.compile(r"DEGREE:\s+(\d+)") - _re_xw = re.compile(r"xw\s?=\s?\[(.+?)\];", re.DOTALL) + re_degree = re.compile(r"DEGREE:\s+(\d+)") + re_xw = re.compile(r"xw\s?=\s?\[(.+?)\];", re.DOTALL) # NOTE: some degrees have multiple quadrature rules with a repeated # header, so we just take the unique ones here degrees = np.unique( - np.fromiter(_re_degree.findall(mfile), dtype=np.int64) + np.fromiter(re_degree.findall(mfile), dtype=np.int64) ) rules = {} - for imatch, match in enumerate(_re_xw.findall(mfile)): + for imatch, match in enumerate(re_xw.findall(mfile)): d = degrees[imatch] assert d == imatch + 1, (d, imatch + 1) From e40277e19d02fbdc16525c72411151f2cf9e55ff Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Thu, 12 Dec 2024 10:07:05 -0600 Subject: [PATCH 08/28] separate quad domain and no quad domain cases; update tests --- modepy/__init__.py | 2 ++ modepy/matrices.py | 45 ++++++++++++++++++++------------- modepy/test/test_matrices.py | 49 +++++++++++++----------------------- 3 files changed, 47 insertions(+), 49 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 8c78884f..49cb9e44 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -38,6 +38,7 @@ multi_vandermonde, nodal_mass_matrix_for_face, nodal_quad_bilinear_form, + nodal_quad_bilinear_form_resampled, nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, resampling_matrix, @@ -167,6 +168,7 @@ "multi_vandermonde", "nodal_mass_matrix_for_face", "nodal_quad_bilinear_form", + "nodal_quad_bilinear_form_resampled", "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", "node_tuples_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 8c672f50..1e311cd4 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -68,6 +68,7 @@ .. autofunction:: spectral_diag_nodal_mass_matrix .. autofunction:: modal_quad_bilinear_form .. autofunction:: nodal_quad_bilinear_form +.. autofunction:: nodal_quad_bilinear_form_resampled Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives @@ -379,8 +380,7 @@ def nodal_quad_bilinear_form( test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], quadrature: Quadrature, - nodes: np.ndarray, - uses_quadrature_domain: bool = False + nodes: np.ndarray ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -391,28 +391,38 @@ def nodal_quad_bilinear_form( 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 point values of the trial - function at either *nodes* or the *quadrature* nodes depending on whether - a quadrature domain is used (as signified by *uses_quadrature_domain*). - - If *uses_quadrature_domain* is set to False, then an interpolation operator - is used to make the operator :math:`N \times N` where :math:`N` is the - number of nodes in *nodes*. Otherwise, the operator has shape :math:`N - \times N_q` where :math:`N_q` is the number of nodes associated with - *quadrature*. Default behavior is to assume a quadrature domain is not used. + function at the *quadrature* nodes. """ if len(test_functions) != nodes.shape[1]: raise ValueError("volume_nodes not unisolvent with trial_functions") vdm_out = vandermonde(trial_functions, nodes) - nodal_operator = la.solve( + return la.solve( vdm_out.T, modal_quad_bilinear_form(quadrature, test_functions)) - if uses_quadrature_domain: - return nodal_operator - else: - return nodal_operator @ resampling_matrix( - trial_functions, quadrature.nodes, nodes) + +def nodal_quad_bilinear_form_resampled( + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + quadrature: Quadrature, + nodes: np.ndarray + ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + + .. math:: + + \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \phi_i(r_k) + \phi_j(r_k)) u_j, + + where :math:`\phi_i` are the Lagrange basis functions obtained from + *test_functions* at *nodes*, :math:`w_k` and :math:`r_k` are the weights + and nodes from *quadrature*, and :math:`u_k` are point values of the trial + function at *nodes*. + """ + return nodal_quad_bilinear_form( + test_functions, trial_functions, quadrature, nodes) @ resampling_matrix( + trial_functions, quadrature.nodes, nodes) def modal_quad_mass_matrix( @@ -462,8 +472,7 @@ def nodal_quad_mass_matrix( test_functions, test_functions, quadrature, - nodes, - uses_quadrature_domain=True) + nodes) def spectral_diag_nodal_mass_matrix( diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index e69d5d23..378b3005 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -49,21 +49,15 @@ def test_nodal_mass_matrix_against_quad( quad_space = mp.space_for_shape(shape, 2*order) quad = mp.quadrature_for_space(quad_space, shape) - if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): - if nodes_on_bdry: - nodes = mp.legendre_gauss_lobatto_tensor_product_nodes( - shape.dim, order, - ) - else: + 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): - if nodes_on_bdry: - nodes = mp.warp_and_blend_nodes(shape.dim, order) - else: + elif isinstance(shape, mp.Simplex): nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes - - else: - raise AssertionError() + else: + raise AssertionError() basis = mp.orthonormal_basis_for_space(space, shape) @@ -136,20 +130,15 @@ def test_bilinear_forms( quad_space = mp.space_for_shape(shape, 2*order) quad = mp.quadrature_for_space(quad_space, shape) - if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): - if nodes_on_bdry: - nodes = mp.legendre_gauss_lobatto_tensor_product_nodes( - shape.dim, order, - ) - else: + 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): - if nodes_on_bdry: - nodes = mp.warp_and_blend_nodes(shape.dim, order) - else: + elif isinstance(shape, mp.Simplex): nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes - else: - raise AssertionError() + else: + raise AssertionError() basis = mp.orthonormal_basis_for_space(space, shape) @@ -160,22 +149,20 @@ def test_bilinear_forms( f = 1 - nodes[ax]**2 fp = -2*nodes[ax] - weak_operator = mp.nodal_quad_bilinear_form( + weak_operator = mp.nodal_quad_bilinear_form_resampled( basis.derivatives(ax), basis.functions, quad, - nodes, - uses_quadrature_domain=False) + nodes) err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) assert err <= 1e-12 else: - quad_mass_mat = mp.nodal_quad_bilinear_form( + quad_mass_mat = mp.nodal_quad_bilinear_form_resampled( basis.functions, basis.functions, quad, - nodes, - uses_quadrature_domain=False) + nodes) vdm_mass_mat = mp.mass_matrix(basis, nodes) err = la.norm(quad_mass_mat - vdm_mass_mat) / la.norm(vdm_mass_mat) From 3881f157c2e62af5441ccfe36e6618f43aeb8bdf Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 11:43:51 -0600 Subject: [PATCH 09/28] minor typo fix --- modepy/matrices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 1e311cd4..b03c7770 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -417,7 +417,7 @@ def nodal_quad_bilinear_form_resampled( where :math:`\phi_i` are the Lagrange basis functions obtained from *test_functions* at *nodes*, :math:`w_k` and :math:`r_k` are the weights - and nodes from *quadrature*, and :math:`u_k` are point values of the trial + and nodes from *quadrature*, and :math:`u_j` are point values of the trial function at *nodes*. """ return nodal_quad_bilinear_form( From 72175ad6c960672a13c78ea27472700a5ad28cb9 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 19:07:10 -0600 Subject: [PATCH 10/28] undo changes to nodal_/modal_quad_mass_matrix --- modepy/matrices.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index b03c7770..38b35c90 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -441,7 +441,12 @@ def modal_quad_mass_matrix( .. versionadded :: 2024.2 """ - return modal_quad_bilinear_form(quadrature, test_functions) + modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) + + for i, test_f in enumerate(test_functions): + modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights + + return modal_mass_matrix def nodal_quad_mass_matrix( @@ -468,11 +473,10 @@ def nodal_quad_mass_matrix( if nodes is None: nodes = quadrature.nodes - return nodal_quad_bilinear_form( - test_functions, - test_functions, - quadrature, - nodes) + vdm = vandermonde(test_functions, nodes) + + return la.solve(vdm.T, + modal_quad_mass_matrix(quadrature, test_functions)) def spectral_diag_nodal_mass_matrix( From c026e7400ec53329d0fa610874cada64319830e3 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 20:35:19 -0600 Subject: [PATCH 11/28] clean things up, make new routine for full bilinear form vs half (?) bilinear form --- modepy/__init__.py | 6 +- modepy/matrices.py | 105 ++++++++++++++++++++++------------- modepy/test/test_matrices.py | 12 ++-- 3 files changed, 77 insertions(+), 46 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 49cb9e44..52bd73c5 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -35,12 +35,13 @@ modal_quad_bilinear_form, modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, + modal_quadrature_operator, multi_vandermonde, nodal_mass_matrix_for_face, nodal_quad_bilinear_form, - nodal_quad_bilinear_form_resampled, nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, + nodal_quadrature_operator, resampling_matrix, spectral_diag_nodal_mass_matrix, vandermonde, @@ -164,13 +165,14 @@ "modal_quad_bilinear_form", "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", + "modal_quadrature_operator", "monomial_basis_for_space", "multi_vandermonde", "nodal_mass_matrix_for_face", "nodal_quad_bilinear_form", - "nodal_quad_bilinear_form_resampled", "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", + "nodal_quadrature_operator", "node_tuples_for_space", "orthonormal_basis_for_space", "quadrature_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 38b35c90..5c8f2b17 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -354,9 +354,9 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) -def modal_quad_bilinear_form( +def modal_quadrature_operator( quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]] + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: r"""Using the *quadrature*, provide a matrix :math:`A` that satisfies: @@ -366,7 +366,8 @@ def modal_quad_bilinear_form( \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`. + :math:`r_j` of *quadrature*, with corresponding weights :math:`w_j`, and + :math:`u_j` is a trial solution evaluated at the *quadrature* nodes. """ modal_operator = np.empty((len(test_functions), len(quadrature.weights))) @@ -376,11 +377,11 @@ def modal_quad_bilinear_form( return modal_operator -def nodal_quad_bilinear_form( - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], +def nodal_quadrature_operator( quadrature: Quadrature, - nodes: np.ndarray + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + nodes: np.ndarray | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -391,38 +392,74 @@ def nodal_quad_bilinear_form( 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 point values of the trial - function at the *quadrature* nodes. + solution at the *quadrature* nodes. """ - if len(test_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with trial_functions") + if nodes is None: + nodes = quadrature.nodes + + if len(proj_functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with interp_functions") - vdm_out = vandermonde(trial_functions, nodes) + vdm = vandermonde(proj_functions, nodes) return la.solve( - vdm_out.T, modal_quad_bilinear_form(quadrature, test_functions)) + vdm.T, modal_quadrature_operator( + quadrature, + test_functions + ) + ) -def nodal_quad_bilinear_form_resampled( - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - quadrature: Quadrature, - nodes: np.ndarray - ) -> np.ndarray: +def modal_quad_bilinear_form( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: .. math:: - \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \phi_i(r_k) + \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \psi_i(r_k) \phi_j(r_k)) u_j, - where :math:`\phi_i` are the Lagrange basis functions obtained from - *test_functions* at *nodes*, :math:`w_k` and :math:`r_k` are the weights - and nodes from *quadrature*, and :math:`u_j` are point values of the trial - function at *nodes*. + where :math:`\psi_i` are from *test_functions*, :math:`\phi_j` are from + *trial_functions*, :math:`r_k` and :math:`w_k` are nodes and weights + from *quadrature*, and :math:`u_j` are modal coefficients of a trial + solution. """ - return nodal_quad_bilinear_form( - test_functions, trial_functions, quadrature, nodes) @ resampling_matrix( - trial_functions, quadrature.nodes, nodes) + return np.einsum( + "qi,qj,q->ij", + vandermonde(test_functions, quadrature.nodes), + vandermonde(trial_functions, quadrature.nodes), + quadrature.weights + ) + + +def nodal_quad_bilinear_form( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + nodes: np.ndarray + ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + + .. math:: + + \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \psi_i(r_k) + \phi_j(r_k)) u_j, + + where :math:`\psi_i` and :math:`\phi_j` are the Lagrange basis functions + obtained from *test_functions* and *trial_functions* at *nodes*, :math:`w_k` + and :math:`r_k` are the weights and nodes from *quadrature*, and :math:`u_j` + are nodal coefficients (point values) of a trial solution at *nodes*. + """ + modal_operator = modal_quad_bilinear_form( + quadrature, test_functions, trial_functions) + + vdm = vandermonde(proj_functions, nodes) + + return la.solve(vdm.T, modal_operator) @ la.inv(vdm) def modal_quad_mass_matrix( @@ -441,12 +478,7 @@ def modal_quad_mass_matrix( .. versionadded :: 2024.2 """ - modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) - - for i, test_f in enumerate(test_functions): - modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights - - return modal_mass_matrix + return modal_quadrature_operator(quadrature, test_functions) def nodal_quad_mass_matrix( @@ -470,13 +502,8 @@ def nodal_quad_mass_matrix( .. versionadded :: 2024.2 """ - if nodes is None: - nodes = quadrature.nodes - - vdm = vandermonde(test_functions, nodes) - - return la.solve(vdm.T, - modal_quad_mass_matrix(quadrature, test_functions)) + return nodal_quadrature_operator( + quadrature, test_functions, test_functions, nodes) def spectral_diag_nodal_mass_matrix( diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 378b3005..88cd41b2 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -114,7 +114,7 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: @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("order", [1, 2, 4]) @pytest.mark.parametrize("nodes_on_bdry", [False, True]) @pytest.mark.parametrize("test_weak_d_dr", [False, True]) def test_bilinear_forms( @@ -149,19 +149,21 @@ def test_bilinear_forms( f = 1 - nodes[ax]**2 fp = -2*nodes[ax] - weak_operator = mp.nodal_quad_bilinear_form_resampled( + weak_operator = mp.nodal_quad_bilinear_form( + quad, basis.derivatives(ax), basis.functions, - quad, + basis.functions, nodes) err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) assert err <= 1e-12 else: - quad_mass_mat = mp.nodal_quad_bilinear_form_resampled( + quad_mass_mat = mp.nodal_quad_bilinear_form( + quad, + basis.functions, basis.functions, basis.functions, - quad, nodes) vdm_mass_mat = mp.mass_matrix(basis, nodes) From fb7fb1fd1bd533b55a9512476eeec566919c023c Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 20:53:39 -0600 Subject: [PATCH 12/28] fix doc failures --- modepy/matrices.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 5c8f2b17..46e6a1ef 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -68,7 +68,8 @@ .. autofunction:: spectral_diag_nodal_mass_matrix .. autofunction:: modal_quad_bilinear_form .. autofunction:: nodal_quad_bilinear_form -.. autofunction:: nodal_quad_bilinear_form_resampled +.. autofunction:: nodal_quadrature_operator +.. autofunction:: modal_quadrature_operator Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives From d3d5ac62c989090fd1576e39ace7fc6d6094f557 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 20:58:32 -0600 Subject: [PATCH 13/28] add order 0 test back to test_bilinear_forms --- modepy/test/test_matrices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 88cd41b2..8b444c0d 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -114,7 +114,7 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None: @pytest.mark.parametrize("shape_cls", [mp.Hypercube, mp.Simplex]) @pytest.mark.parametrize("dim", [1, 2, 3]) -@pytest.mark.parametrize("order", [1, 2, 4]) +@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( From c063d234325d04823773c3ea614ee7f2e0bb8586 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 13 Dec 2024 21:03:34 -0600 Subject: [PATCH 14/28] typo fix, add an additional assert in nodal_bilinear_forms --- modepy/matrices.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 46e6a1ef..061f263b 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -399,7 +399,7 @@ def nodal_quadrature_operator( nodes = quadrature.nodes if len(proj_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with interp_functions") + raise ValueError("volume_nodes not unisolvent with proj_functions") vdm = vandermonde(proj_functions, nodes) @@ -455,6 +455,9 @@ def nodal_quad_bilinear_form( and :math:`r_k` are the weights and nodes from *quadrature*, and :math:`u_j` are nodal coefficients (point values) of a trial solution at *nodes*. """ + if len(test_functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with proj_functions") + modal_operator = modal_quad_bilinear_form( quadrature, test_functions, trial_functions) From d3f9e365a200f760ae2c45a73a7dfcfbba96df4e Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Sat, 14 Dec 2024 11:22:00 -0600 Subject: [PATCH 15/28] *_quadrature_operator -> *_quad_operator --- modepy/__init__.py | 8 ++++---- modepy/matrices.py | 20 ++++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 52bd73c5..ce9e3fa2 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -35,13 +35,13 @@ modal_quad_bilinear_form, modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, - modal_quadrature_operator, + modal_quad_operator, multi_vandermonde, nodal_mass_matrix_for_face, nodal_quad_bilinear_form, nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, - nodal_quadrature_operator, + nodal_quad_operator, resampling_matrix, spectral_diag_nodal_mass_matrix, vandermonde, @@ -165,14 +165,14 @@ "modal_quad_bilinear_form", "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", - "modal_quadrature_operator", + "modal_quad_operator", "monomial_basis_for_space", "multi_vandermonde", "nodal_mass_matrix_for_face", "nodal_quad_bilinear_form", "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", - "nodal_quadrature_operator", + "nodal_quad_operator", "node_tuples_for_space", "orthonormal_basis_for_space", "quadrature_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 061f263b..0214f257 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -68,8 +68,8 @@ .. autofunction:: spectral_diag_nodal_mass_matrix .. autofunction:: modal_quad_bilinear_form .. autofunction:: nodal_quad_bilinear_form -.. autofunction:: nodal_quadrature_operator -.. autofunction:: modal_quadrature_operator +.. autofunction:: nodal_quad_operator +.. autofunction:: modal_quad_operator Differentiation is also convenient to express by using :math:`V^{-1}` to obtain modal values and then using a Vandermonde matrix for the derivatives @@ -355,7 +355,7 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) -def modal_quadrature_operator( +def modal_quad_operator( quadrature: Quadrature, test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: @@ -378,7 +378,7 @@ def modal_quadrature_operator( return modal_operator -def nodal_quadrature_operator( +def nodal_quad_operator( quadrature: Quadrature, test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], @@ -404,7 +404,7 @@ def nodal_quadrature_operator( vdm = vandermonde(proj_functions, nodes) return la.solve( - vdm.T, modal_quadrature_operator( + vdm.T, modal_quad_operator( quadrature, test_functions ) @@ -416,12 +416,12 @@ def modal_quad_bilinear_form( test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + r"""Using *quadrature*, provide a matrix :math:`A` defined as: .. math:: - \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \psi_i(r_k) - \phi_j(r_k)) u_j, + \displaystyle A_{i,j} = \sum_k w_k \psi_i(r_k) + \phi_j(r_k), where :math:`\psi_i` are from *test_functions*, :math:`\phi_j` are from *trial_functions*, :math:`r_k` and :math:`w_k` are nodes and weights @@ -482,7 +482,7 @@ def modal_quad_mass_matrix( .. versionadded :: 2024.2 """ - return modal_quadrature_operator(quadrature, test_functions) + return modal_quad_operator(quadrature, test_functions) def nodal_quad_mass_matrix( @@ -506,7 +506,7 @@ def nodal_quad_mass_matrix( .. versionadded :: 2024.2 """ - return nodal_quadrature_operator( + return nodal_quad_operator( quadrature, test_functions, test_functions, nodes) From 6adb3e78a599ce470496f941a64627cf22993c8c Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Sat, 14 Dec 2024 13:41:19 -0600 Subject: [PATCH 16/28] clarify docs --- modepy/matrices.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 0214f257..fccdd9a5 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -420,13 +420,17 @@ def modal_quad_bilinear_form( .. math:: - \displaystyle A_{i,j} = \sum_k w_k \psi_i(r_k) - \phi_j(r_k), + \displaystyle A_{ij} = \sum_k \psi_i(r_k) \phi_j(r_k) w_k, where :math:`\psi_i` are from *test_functions*, :math:`\phi_j` are from *trial_functions*, :math:`r_k` and :math:`w_k` are nodes and weights - from *quadrature*, and :math:`u_j` are modal coefficients of a trial - solution. + from *quadrature*. The matrix :math:`A` satisfies + + .. math:: + + \displaystyle (u, psi_i)_A = \sum_{j} A_{ij} u_j, + + where :math:`u_i` are modal coefficients of a trial solution. """ return np.einsum( "qi,qj,q->ij", @@ -443,17 +447,22 @@ def nodal_quad_bilinear_form( proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], nodes: np.ndarray ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: + r"""Using *quadrature*, provide a matrix :math:`A` defined as: + + .. math:: + + \displaystyle A_{ij} = \sum_k \psi_i(r_k) \phi_j(r_k) w_k, + + 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 .. math:: - \displaystyle (A \boldsymbol u)_i = \sum_k,j (w_k \psi_i(r_k) - \phi_j(r_k)) u_j, + \displaystyle (u, psi_i)_A = \sum_{j} A_{ij} u_j, - where :math:`\psi_i` and :math:`\phi_j` are the Lagrange basis functions - obtained from *test_functions* and *trial_functions* at *nodes*, :math:`w_k` - and :math:`r_k` are the weights and nodes from *quadrature*, and :math:`u_j` - are nodal coefficients (point values) of a trial solution at *nodes*. + where :math:`u_i` are nodal coefficients of a trial solution. """ if len(test_functions) != nodes.shape[1]: raise ValueError("volume_nodes not unisolvent with proj_functions") From f914a0a91c78c522b450356142ad7e8826836374 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Sat, 14 Dec 2024 13:46:01 -0600 Subject: [PATCH 17/28] more doc clarification --- modepy/matrices.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index fccdd9a5..07a49515 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -359,7 +359,7 @@ def modal_quad_operator( quadrature: Quadrature, test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: - r"""Using the *quadrature*, provide a matrix :math:`A` that + r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: .. math:: @@ -392,7 +392,7 @@ def nodal_quad_operator( 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 point values of the trial + and nodes from *quadrature*, and :math:`u_j` are point values of a trial solution at the *quadrature* nodes. """ if nodes is None: @@ -428,7 +428,7 @@ def modal_quad_bilinear_form( .. math:: - \displaystyle (u, psi_i)_A = \sum_{j} A_{ij} u_j, + \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, where :math:`u_i` are modal coefficients of a trial solution. """ @@ -460,7 +460,7 @@ def nodal_quad_bilinear_form( .. math:: - \displaystyle (u, psi_i)_A = \sum_{j} A_{ij} u_j, + \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, where :math:`u_i` are nodal coefficients of a trial solution. """ @@ -479,7 +479,7 @@ 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 + r"""Using *quadrature*, provide a matrix :math:`M` that satisfies: .. math:: @@ -499,7 +499,7 @@ def nodal_quad_mass_matrix( test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], nodes: np.ndarray | None = None, ) -> np.ndarray: - r"""Using the *quadrature*, provide a matrix :math:`M` that + r"""Using *quadrature*, provide a matrix :math:`M` that satisfies: .. math:: From 7b82ab5b568b26d8b5eac910887ea76c489865b9 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Mon, 23 Dec 2024 15:35:34 -0600 Subject: [PATCH 18/28] deprecate old functions, delete tests, beef up nodal bilinear form --- modepy/__init__.py | 2 - modepy/matrices.py | 237 +++++++++++++---------------------- modepy/test/test_matrices.py | 58 ++------- 3 files changed, 99 insertions(+), 198 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index ce9e3fa2..194a29d0 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -32,7 +32,6 @@ inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face, - modal_quad_bilinear_form, modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, modal_quad_operator, @@ -162,7 +161,6 @@ "legendre_gauss_tensor_product_nodes", "mass_matrix", "modal_mass_matrix_for_face", - "modal_quad_bilinear_form", "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", "modal_quad_operator", diff --git a/modepy/matrices.py b/modepy/matrices.py index 07a49515..555366e0 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -82,8 +82,11 @@ """ +BasisFunctionType = Callable[[np.ndarray], np.ndarray] + + def vandermonde( - functions: Sequence[Callable[[np.ndarray], np.ndarray]], + functions: Sequence[BasisFunctionType], nodes: np.ndarray ) -> np.ndarray: """Return a (generalized) Vandermonde matrix. @@ -167,7 +170,7 @@ def multi_vandermonde( def resampling_matrix( - basis: Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Sequence[BasisFunctionType], new_nodes: np.ndarray, old_nodes: np.ndarray, least_squares_ok: bool = False @@ -222,7 +225,7 @@ def is_square(m): def differentiation_matrices( - basis: Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Sequence[BasisFunctionType], grad_basis: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]], nodes: np.ndarray, from_nodes: np.ndarray | None = None @@ -302,7 +305,7 @@ def diff_matrix_permutation( def inverse_mass_matrix( - basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Basis | Sequence[BasisFunctionType], nodes: np.ndarray ) -> np.ndarray: """Return a matrix :math:`A=M^{-1}`, which is the inverse of the one returned @@ -321,7 +324,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[BasisFunctionType] = ( basis.functions) else: basis_functions = basis @@ -337,7 +340,7 @@ def inverse_mass_matrix( def mass_matrix( - basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]], + basis: Basis | Sequence[BasisFunctionType], nodes: np.ndarray ) -> np.ndarray: r"""Return a mass matrix :math:`M`, which obeys @@ -357,7 +360,7 @@ def mass_matrix( def modal_quad_operator( quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + test_functions: Sequence[BasisFunctionType], ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -380,9 +383,10 @@ def modal_quad_operator( def nodal_quad_operator( quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - nodes: np.ndarray | None = None + test_functions: Sequence[BasisFunctionType], + proj_functions: Sequence[BasisFunctionType], + nodes: np.ndarray | None = None, + mapped_nodes: np.ndarray | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -411,41 +415,15 @@ def nodal_quad_operator( ) -def modal_quad_bilinear_form( - quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` defined as: - - .. math:: - - \displaystyle A_{ij} = \sum_k \psi_i(r_k) \phi_j(r_k) w_k, - - where :math:`\psi_i` are from *test_functions*, :math:`\phi_j` are from - *trial_functions*, :math:`r_k` and :math:`w_k` are nodes and weights - from *quadrature*. The matrix :math:`A` satisfies - - .. math:: - - \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, - - where :math:`u_i` are modal coefficients of a trial solution. - """ - return np.einsum( - "qi,qj,q->ij", - vandermonde(test_functions, quadrature.nodes), - vandermonde(trial_functions, quadrature.nodes), - quadrature.weights - ) - - def nodal_quad_bilinear_form( quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - proj_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - nodes: np.ndarray + test_basis: Basis, + trial_basis: Basis, + input_nodes: np.ndarray, + output_nodes: np.ndarray | None = None, + mapping_function: Callable | None = None, + test_derivative_ax: int | None = None, + trial_derivative_ax: int | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` defined as: @@ -464,59 +442,35 @@ def nodal_quad_bilinear_form( where :math:`u_i` are nodal coefficients of a trial solution. """ - if len(test_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with proj_functions") - - modal_operator = modal_quad_bilinear_form( - quadrature, test_functions, trial_functions) - - vdm = vandermonde(proj_functions, nodes) - - return la.solve(vdm.T, modal_operator) @ la.inv(vdm) - - -def modal_quad_mass_matrix( - quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`M` that - satisfies: - - .. math:: - - \displaystyle (M \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`. - - .. versionadded :: 2024.2 - """ - return modal_quad_operator(quadrature, test_functions) + if output_nodes is None: + output_nodes = input_nodes + test_functions = ( + test_basis.derivatives(test_derivative_ax) + if test_derivative_ax is not None else test_basis.functions + ) -def nodal_quad_mass_matrix( - quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - nodes: np.ndarray | None = None, - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`M` that - satisfies: - - .. math:: + trial_functions = ( + trial_basis.derivatives(trial_derivative_ax) + if trial_derivative_ax is not None else trial_basis.functions + ) - \displaystyle (M \boldsymbol u)_i = \sum_j w_j \phi_i(r_j) u_j, + mapped_nodes = ( + mapping_function(quadrature.nodes) + if mapping_function is not None else quadrature.nodes + ) - 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*. + modal_operator = np.einsum( + "qi,qj,q->ij", + vandermonde(test_functions, mapped_nodes), + vandermonde(trial_functions, quadrature.nodes), + quadrature.weights + ) - If *nodes* is not provided, use *quadrature*'s nodes. + input_vdm = vandermonde(trial_basis.functions, input_nodes) + output_vdm = vandermonde(test_basis.functions, output_nodes) - .. versionadded :: 2024.2 - """ - return nodal_quad_operator( - quadrature, test_functions, test_functions, nodes) + return la.solve(output_vdm.T, modal_operator) @ la.inv(input_vdm) def spectral_diag_nodal_mass_matrix( @@ -540,25 +494,38 @@ 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]] - ) -> np.ndarray: - r"""Using the quadrature *face_quad*, provide a matrix :math:`M_f` that - satisfies: +# {{{ deprecated remove in 2025-ish - .. math:: +def modal_quad_mass_matrix( + quadrature: Quadrature, + test_functions: Sequence[BasisFunctionType], + ) -> np.ndarray: + from warnings import warn + warn("`modal_quad_mass_matrix` is deprecated and will stop working in " + "2025. Please use `modal_quad_operator` instead.", stacklevel=1) + return modal_quad_operator(quadrature, test_functions) - \displaystyle {(M_f)}_{ij} \approx \int_F \phi_i(r) \psi_j(r) dr, - where :math:`\phi_i` are the (volume) basis functions - *test_functions*, and :math:`\psi_j` are the (surface) - *trial_functions*. +def nodal_quad_mass_matrix( + quadrature: Quadrature, + test_functions: Sequence[BasisFunctionType], + nodes: np.ndarray | None = None, + ) -> np.ndarray: + from warnings import warn + warn("`nodal_quad_mass_matrix` is deprecated and will stop working in " + "2025. Please use `nodal_quad_operator` instead.", stacklevel=1) + return nodal_quad_operator( + quadrature, test_functions, test_functions, nodes) - .. versionadded:: 2020.3 - """ +def modal_mass_matrix_for_face( + face: Face, face_quad: Quadrature, + trial_functions: Sequence[BasisFunctionType], + test_functions: Sequence[BasisFunctionType] + ) -> np.ndarray: + from warnings import warn + warn("`modal_mass_matrix_for_face` is deprecated and will stop working in " + "2025. Please use `modal_quad_bilinear_form` instead.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) result = np.empty((len(test_functions), len(trial_functions))) @@ -573,23 +540,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[BasisFunctionType], + test_functions: Sequence[BasisFunctionType], + 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 " + "2025. Please use `nodal_quad_bilinear_form` instead.", stacklevel=1) face_vdm = vandermonde(trial_functions, face_nodes) vol_vdm = vandermonde(test_functions, volume_nodes) @@ -608,21 +566,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[BasisFunctionType], ) -> 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. Please use `modal_quad_operator` instead.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) vol_modal_mass_matrix = np.empty((len(test_functions), len(face_quad.weights))) @@ -635,23 +583,12 @@ 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[BasisFunctionType], 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. Please use `nodal_quad_operator` instead.", stacklevel=1) if len(test_functions) != volume_nodes.shape[1]: raise ValueError("volume_nodes not unisolvent with test_functions") @@ -660,4 +597,6 @@ def nodal_quad_mass_matrix_for_face( return la.solve(vol_vdm.T, modal_quad_mass_matrix_for_face(face, face_quad, test_functions)) +# }}} + # vim: foldmethod=marker diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 8b444c0d..9ce6d78a 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -33,43 +33,6 @@ import modepy as mp -@pytest.mark.parametrize("shape", [ - mp.Simplex(1), - mp.Simplex(2), - mp.Simplex(3), - mp.Hypercube(2), - mp.Hypercube(3), - ]) -@pytest.mark.parametrize("order", [0, 1, 2, 4]) -@pytest.mark.parametrize("nodes_on_bdry", [False, True]) -def test_nodal_mass_matrix_against_quad( - shape: mp.Shape, order: int, nodes_on_bdry: bool, - ) -> None: - 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) - - quad_mass_mat = mp.nodal_quad_mass_matrix(quad, basis.functions, nodes) - vdm_mass_mat = mp.mass_matrix(basis, nodes) - - nodes_to_quad = mp.resampling_matrix(basis.functions, quad.nodes, nodes) - - err = la.norm(quad_mass_mat@nodes_to_quad - vdm_mass_mat)/la.norm(vdm_mass_mat) - assert err < 1e-14 - - @pytest.mark.parametrize("shape", [ mp.Hypercube(1), mp.Hypercube(2), @@ -150,21 +113,22 @@ def test_bilinear_forms( fp = -2*nodes[ax] weak_operator = mp.nodal_quad_bilinear_form( - quad, - basis.derivatives(ax), - basis.functions, - basis.functions, - nodes) + quadrature=quad, + test_basis=basis, + trial_basis=basis, + input_nodes=nodes, + test_derivative_ax=ax + ) err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) assert err <= 1e-12 else: quad_mass_mat = mp.nodal_quad_bilinear_form( - quad, - basis.functions, - basis.functions, - basis.functions, - nodes) + quadrature=quad, + test_basis=basis, + trial_basis=basis, + 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) From e5c496cc24975db4cc11359499734a7393d1ca39 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Mon, 23 Dec 2024 15:44:08 -0600 Subject: [PATCH 19/28] restore nodal/modal quad mass (keep deprecated) --- modepy/matrices.py | 68 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 20 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 555366e0..b145562b 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -415,6 +415,25 @@ def nodal_quad_operator( ) +def modal_quad_bilinear_form( + quadrature: Quadrature, + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], + mapping_function: Callable[[np.ndarray], np.ndarray] | None = None, + ) -> np.ndarray: + mapped_nodes = ( + mapping_function(quadrature.nodes) + if mapping_function is not None else quadrature.nodes + ) + + return np.einsum( + "qi,qj,q->ij", + vandermonde(test_functions, mapped_nodes), + vandermonde(trial_functions, quadrature.nodes), + quadrature.weights + ) + + def nodal_quad_bilinear_form( quadrature: Quadrature, test_basis: Basis, @@ -455,16 +474,11 @@ def nodal_quad_bilinear_form( if trial_derivative_ax is not None else trial_basis.functions ) - mapped_nodes = ( - mapping_function(quadrature.nodes) - if mapping_function 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 + modal_operator = modal_quad_bilinear_form( + quadrature, + test_functions, + trial_functions, + mapping_function ) input_vdm = vandermonde(trial_basis.functions, input_nodes) @@ -498,24 +512,37 @@ def spectral_diag_nodal_mass_matrix( def modal_quad_mass_matrix( quadrature: Quadrature, - test_functions: Sequence[BasisFunctionType], + test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], ) -> np.ndarray: from warnings import warn warn("`modal_quad_mass_matrix` is deprecated and will stop working in " - "2025. Please use `modal_quad_operator` instead.", stacklevel=1) - return modal_quad_operator(quadrature, test_functions) + "2025.", stacklevel=0) + modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) + + for i, test_f in enumerate(test_functions): + modal_mass_matrix[i] = test_f(quadrature.nodes) * quadrature.weights + + return modal_mass_matrix def nodal_quad_mass_matrix( quadrature: Quadrature, - test_functions: Sequence[BasisFunctionType], + 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 " - "2025. Please use `nodal_quad_operator` instead.", stacklevel=1) - return nodal_quad_operator( - quadrature, test_functions, test_functions, nodes) + "2025. Consider switching to `nodal_quad_bilinear_form`", stacklevel=0) + 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( @@ -525,7 +552,7 @@ def modal_mass_matrix_for_face( ) -> np.ndarray: from warnings import warn warn("`modal_mass_matrix_for_face` is deprecated and will stop working in " - "2025. Please use `modal_quad_bilinear_form` instead.", stacklevel=1) + "2025.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) result = np.empty((len(test_functions), len(trial_functions))) @@ -570,7 +597,7 @@ def modal_quad_mass_matrix_for_face( ) -> np.ndarray: from warnings import warn warn("`modal_quad_mass_matrix_for_face` is deprecated and will stop working " - "in 2025. Please use `modal_quad_operator` instead.", stacklevel=1) + "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))) @@ -588,7 +615,8 @@ def nodal_quad_mass_matrix_for_face( ) -> np.ndarray: from warnings import warn warn("`nodal_quad_mass_matrix_for_face` is deprecated and will stop working " - "in 2025. Please use `nodal_quad_operator` instead.", stacklevel=1) + "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") From 6b2f0537e10ce9a5b2284073eae693edd0b06df9 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Mon, 23 Dec 2024 15:54:09 -0600 Subject: [PATCH 20/28] update docs --- modepy/matrices.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/modepy/matrices.py b/modepy/matrices.py index b145562b..caf101f2 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -421,6 +421,21 @@ def modal_quad_bilinear_form( trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], mapping_function: Callable[[np.ndarray], np.ndarray] | None = None, ) -> np.ndarray: + r"""Using *quadrature*, provide a matrix :math:`A` defined as: + + .. math:: + + \displaystyle A_{ij} = \sum_k \phi_i(r_k) \psi_j(r_k) w_k, + + where :math:`phi_i` are in *test_functions*, :math:`\psi_j` are in + *trial_functions*, and :math:`r_k` and :math:`w_k` are nodes and weights + from *quadrature*. + + An optional *mapping_function* can be supplied to map the arguments + evaluated by *test_functions*. This would be useful, e.g., to evaluate a + face mass matrix in the context of DG-FEM. Otherwise, the input to + *test_functions* are the *quadrature* nodes. + """ mapped_nodes = ( mapping_function(quadrature.nodes) if mapping_function is not None else quadrature.nodes @@ -460,6 +475,11 @@ def nodal_quad_bilinear_form( \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, where :math:`u_i` are nodal coefficients of a trial solution. + + Optional arguments allow this function to work in cases where the input and + output discretizations are not necessarily the same. This is the case, for + example, in DG-FEM when evaluating a face mass matrix or when there is a + base and quadrature discretization. """ if output_nodes is None: output_nodes = input_nodes From e1b5ae1d8173b4c5702401f544155e669b3e01ae Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Tue, 31 Dec 2024 00:05:17 -0600 Subject: [PATCH 21/28] move facemass to use general routines --- modepy/matrices.py | 34 +++++++++++++++++---------- modepy/test/test_tools.py | 48 +++++++++++++++++++++++++++------------ 2 files changed, 56 insertions(+), 26 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index caf101f2..336722c1 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -383,10 +383,10 @@ def modal_quad_operator( def nodal_quad_operator( quadrature: Quadrature, - test_functions: Sequence[BasisFunctionType], - proj_functions: Sequence[BasisFunctionType], + test_basis: Basis, + test_derivative_ax: int | None = None, nodes: np.ndarray | None = None, - mapped_nodes: np.ndarray | None = None + mapping_function: Callable[[np.ndarray], np.ndarray] | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -402,18 +402,28 @@ def nodal_quad_operator( if nodes is None: nodes = quadrature.nodes - if len(proj_functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with proj_functions") + if len(test_basis.functions) != nodes.shape[1]: + raise ValueError("volume_nodes not unisolvent with test functions") + + test_functions = ( + test_basis.derivatives(test_derivative_ax) + if test_derivative_ax is not None else test_basis.functions + ) - vdm = vandermonde(proj_functions, nodes) + vdm = vandermonde(test_basis.functions, nodes) - return la.solve( - vdm.T, modal_quad_operator( - quadrature, - test_functions - ) + mapped_nodes = ( + mapping_function(quadrature.nodes) + if mapping_function is not None else quadrature.nodes ) + modal_operator = np.array([ + test_f(mapped_nodes) * quadrature.weights + for test_f in test_functions + ]) + + return la.solve(vdm.T, modal_operator) + def modal_quad_bilinear_form( quadrature: Quadrature, @@ -504,7 +514,7 @@ def nodal_quad_bilinear_form( input_vdm = vandermonde(trial_basis.functions, input_nodes) output_vdm = vandermonde(test_basis.functions, output_nodes) - return la.solve(output_vdm.T, modal_operator) @ la.inv(input_vdm) + return la.solve(output_vdm.T, modal_operator @ la.inv(input_vdm)) def spectral_diag_nodal_mass_matrix( diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py index ef575661..df49ffb7 100644 --- a/modepy/test/test_tools.py +++ b/modepy/test/test_tools.py @@ -295,22 +295,31 @@ 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): +def test_modal_quad_bilinear_form_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 + from modepy.matrices import modal_quad_bilinear_form 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) + fmm = modal_quad_bilinear_form( + quadrature=face_quad, + test_functions=vol_basis.functions, + trial_functions=face_basis.functions, + mapping_function=face.map_to_volume + ) + + fmm2 = modal_quad_bilinear_form( + quadrature=face_quad2, + test_functions=vol_basis.functions, + trial_functions=face_basis.functions, + mapping_function=face.map_to_volume + ) error = la.norm(fmm - fmm2, np.inf) / la.norm(fmm2, np.inf) logger.info("fmm error: %.5e", error) @@ -324,7 +333,7 @@ def test_modal_mass_matrix_for_face(dims, shape_cls, order=3): @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_quad_bilinear_form_for_face(dims, shape_cls, order=3): vol_shape = shape_cls(dims) vol_space = mp.space_for_shape(vol_shape, order) @@ -332,8 +341,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_quad_bilinear_form, + nodal_quad_operator, ) for face in mp.faces_for_shape(vol_shape): face_space = mp.space_for_shape(face, order) @@ -342,11 +351,22 @@ 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_quad_bilinear_form( + quadrature=face_quad, + test_basis=volume_basis, + trial_basis=face_basis, + input_nodes=face_nodes, + output_nodes=volume_nodes, + mapping_function=face.map_to_volume + ) + + fmm2 = nodal_quad_operator( + quadrature=face_quad2, + test_basis=volume_basis, + nodes=volume_nodes, + mapping_function=face.map_to_volume + ) for f_face in face_basis.functions: fval_nodal = f_face(face_nodes) From 6157efc41738f40555eccefaf26f1ed6ab85b1e0 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Wed, 29 Jan 2025 10:17:03 -0600 Subject: [PATCH 22/28] update docs, remove modal routines, deprecate mass, face mass routines --- modepy/__init__.py | 14 ++-- modepy/matrices.py | 121 +++++++++++------------------------ modepy/test/test_matrices.py | 4 +- modepy/test/test_tools.py | 48 ++------------ 4 files changed, 47 insertions(+), 140 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 194a29d0..febbe22c 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -32,15 +32,12 @@ inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face, - modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, - modal_quad_operator, multi_vandermonde, nodal_mass_matrix_for_face, - nodal_quad_bilinear_form, - nodal_quad_mass_matrix, nodal_quad_mass_matrix_for_face, - nodal_quad_operator, + nodal_quadrature_bilinear_form_matrix, + nodal_quadrature_matrix, resampling_matrix, spectral_diag_nodal_mass_matrix, vandermonde, @@ -161,16 +158,13 @@ "legendre_gauss_tensor_product_nodes", "mass_matrix", "modal_mass_matrix_for_face", - "modal_quad_mass_matrix", "modal_quad_mass_matrix_for_face", - "modal_quad_operator", "monomial_basis_for_space", "multi_vandermonde", "nodal_mass_matrix_for_face", - "nodal_quad_bilinear_form", - "nodal_quad_mass_matrix", "nodal_quad_mass_matrix_for_face", - "nodal_quad_operator", + "nodal_quadrature_bilinear_form_matrix", + "nodal_quadrature_matrix", "node_tuples_for_space", "orthonormal_basis_for_space", "quadrature_for_space", diff --git a/modepy/matrices.py b/modepy/matrices.py index 336722c1..83dee9f3 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -61,15 +61,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:: modal_quad_bilinear_form -.. autofunction:: nodal_quad_bilinear_form -.. autofunction:: nodal_quad_operator -.. autofunction:: modal_quad_operator +.. autofunction:: nodal_quadrature_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 @@ -358,30 +352,7 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) -def modal_quad_operator( - quadrature: Quadrature, - test_functions: Sequence[BasisFunctionType], - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` that - satisfies: - - .. math:: - - \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`, and - :math:`u_j` is a trial solution evaluated at the *quadrature* nodes. - """ - modal_operator = np.empty((len(test_functions), len(quadrature.weights))) - - for i, test_f in enumerate(test_functions): - modal_operator[i] = test_f(quadrature.nodes) * quadrature.weights - - return modal_operator - - -def nodal_quad_operator( +def nodal_quadrature_matrix( quadrature: Quadrature, test_basis: Basis, test_derivative_ax: int | None = None, @@ -396,8 +367,8 @@ def nodal_quad_operator( 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 point values of a trial - solution at the *quadrature* nodes. + and nodes from *quadrature*, and :math:`u_j` are trial solution point values + at *quadrature* nodes. """ if nodes is None: nodes = quadrature.nodes @@ -425,47 +396,13 @@ def nodal_quad_operator( return la.solve(vdm.T, modal_operator) -def modal_quad_bilinear_form( - quadrature: Quadrature, - test_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]], - mapping_function: Callable[[np.ndarray], np.ndarray] | None = None, - ) -> np.ndarray: - r"""Using *quadrature*, provide a matrix :math:`A` defined as: - - .. math:: - - \displaystyle A_{ij} = \sum_k \phi_i(r_k) \psi_j(r_k) w_k, - - where :math:`phi_i` are in *test_functions*, :math:`\psi_j` are in - *trial_functions*, and :math:`r_k` and :math:`w_k` are nodes and weights - from *quadrature*. - - An optional *mapping_function* can be supplied to map the arguments - evaluated by *test_functions*. This would be useful, e.g., to evaluate a - face mass matrix in the context of DG-FEM. Otherwise, the input to - *test_functions* are the *quadrature* nodes. - """ - mapped_nodes = ( - mapping_function(quadrature.nodes) - if mapping_function is not None else quadrature.nodes - ) - - return np.einsum( - "qi,qj,q->ij", - vandermonde(test_functions, mapped_nodes), - vandermonde(trial_functions, quadrature.nodes), - quadrature.weights - ) - - -def nodal_quad_bilinear_form( +def nodal_quadrature_bilinear_form_matrix( quadrature: Quadrature, test_basis: Basis, trial_basis: Basis, input_nodes: np.ndarray, output_nodes: np.ndarray | None = None, - mapping_function: Callable | None = None, + mapping_function: Callable[[np.ndarray], np.ndarray] | None = None, test_derivative_ax: int | None = None, trial_derivative_ax: int | None = None ) -> np.ndarray: @@ -484,12 +421,20 @@ def nodal_quad_bilinear_form( \displaystyle (u, \psi_i)_A = \sum_{j} A_{ij} u_j, - where :math:`u_i` are nodal coefficients of a trial solution. + where :math:`u_i` are nodal coefficients. - Optional arguments allow this function to work in cases where the input and - output discretizations are not necessarily the same. This is the case, for - example, in DG-FEM when evaluating a face mass matrix or when there is a - base and quadrature discretization. + An optional *mapping_function* can be supplied to map the *quadrature* nodes + evaluated by *test_basis* before constructing the matrix. + + An optional set of *output_nodes* can be supplied if the result should live + on a domain different from the domain of *input_nodes*. + + Test and trial function derivatives can be requested by supplying either + *test_derivative_ax* or *trial_derivative_ax* as an integer corresponding to + an enumerated spatial derivative ax. For example, *test_derivative_ax = 0* + would use *test_basis* differentiated with respect to :math:`x` while + *trial_derivative_ax = 2* would use *trial_basis* differentiated with + respect to :math:`z`. """ if output_nodes is None: output_nodes = input_nodes @@ -504,11 +449,16 @@ def nodal_quad_bilinear_form( if trial_derivative_ax is not None else trial_basis.functions ) - modal_operator = modal_quad_bilinear_form( - quadrature, - test_functions, - trial_functions, - mapping_function + mapped_nodes = ( + mapping_function(quadrature.nodes) + if mapping_function 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(trial_basis.functions, input_nodes) @@ -546,7 +496,7 @@ def modal_quad_mass_matrix( ) -> np.ndarray: from warnings import warn warn("`modal_quad_mass_matrix` is deprecated and will stop working in " - "2025.", stacklevel=0) + "2026.", stacklevel=1) modal_mass_matrix = np.empty((len(test_functions), len(quadrature.weights))) for i, test_f in enumerate(test_functions): @@ -562,7 +512,7 @@ def nodal_quad_mass_matrix( ) -> np.ndarray: from warnings import warn warn("`nodal_quad_mass_matrix` is deprecated and will stop working in " - "2025. Consider switching to `nodal_quad_bilinear_form`", stacklevel=0) + "2026. Consider switching to `nodal_quad_bilinear_form`", stacklevel=1) if nodes is None: nodes = quadrature.nodes @@ -582,7 +532,7 @@ def modal_mass_matrix_for_face( ) -> np.ndarray: from warnings import warn warn("`modal_mass_matrix_for_face` is deprecated and will stop working in " - "2025.", stacklevel=1) + "2026.", stacklevel=1) mapped_nodes = face.map_to_volume(face_quad.nodes) result = np.empty((len(test_functions), len(trial_functions))) @@ -604,7 +554,7 @@ def nodal_mass_matrix_for_face( ) -> np.ndarray: from warnings import warn warn("`nodal_mass_matrix_for_face` is deprecated and will stop working in " - "2025. Please use `nodal_quad_bilinear_form` instead.", stacklevel=1) + "2026. Please use `nodal_quad_bilinear_form` instead.", stacklevel=1) face_vdm = vandermonde(trial_functions, face_nodes) vol_vdm = vandermonde(test_functions, volume_nodes) @@ -653,7 +603,8 @@ def nodal_quad_mass_matrix_for_face( 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)) # }}} diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 9ce6d78a..ca130925 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -112,7 +112,7 @@ def test_bilinear_forms( f = 1 - nodes[ax]**2 fp = -2*nodes[ax] - weak_operator = mp.nodal_quad_bilinear_form( + weak_operator = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, test_basis=basis, trial_basis=basis, @@ -123,7 +123,7 @@ def test_bilinear_forms( err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) assert err <= 1e-12 else: - quad_mass_mat = mp.nodal_quad_bilinear_form( + quad_mass_mat = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, test_basis=basis, trial_basis=basis, diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py index df49ffb7..8dc8d437 100644 --- a/modepy/test/test_tools.py +++ b/modepy/test/test_tools.py @@ -295,45 +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_quad_bilinear_form_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_quad_bilinear_form - 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_quad_bilinear_form( - quadrature=face_quad, - test_functions=vol_basis.functions, - trial_functions=face_basis.functions, - mapping_function=face.map_to_volume - ) - - fmm2 = modal_quad_bilinear_form( - quadrature=face_quad2, - test_functions=vol_basis.functions, - trial_functions=face_basis.functions, - mapping_function=face.map_to_volume - ) - - 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_quad_bilinear_form_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) @@ -341,8 +303,8 @@ def test_nodal_quad_bilinear_form_for_face(dims, shape_cls, order=3): volume_basis = mp.basis_for_space(vol_space, vol_shape) from modepy.matrices import ( - nodal_quad_bilinear_form, - nodal_quad_operator, + nodal_quadrature_bilinear_form_matrix, + nodal_quadrature_matrix, ) for face in mp.faces_for_shape(vol_shape): face_space = mp.space_for_shape(face, order) @@ -352,7 +314,7 @@ def test_nodal_quad_bilinear_form_for_face(dims, shape_cls, order=3): face_quad2 = mp.quadrature_for_space( mp.space_for_shape(face, 2*order+2), face) - fmm = nodal_quad_bilinear_form( + fmm = nodal_quadrature_bilinear_form_matrix( quadrature=face_quad, test_basis=volume_basis, trial_basis=face_basis, @@ -361,7 +323,7 @@ def test_nodal_quad_bilinear_form_for_face(dims, shape_cls, order=3): mapping_function=face.map_to_volume ) - fmm2 = nodal_quad_operator( + fmm2 = nodal_quadrature_matrix( quadrature=face_quad2, test_basis=volume_basis, nodes=volume_nodes, From e117c24c7e312dac09f0b03d783a1f7118b5ed52 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Wed, 29 Jan 2025 10:19:40 -0600 Subject: [PATCH 23/28] fix typos, add unisolvency check --- modepy/matrices.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 83dee9f3..e9a7d153 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -374,7 +374,7 @@ def nodal_quadrature_matrix( nodes = quadrature.nodes if len(test_basis.functions) != nodes.shape[1]: - raise ValueError("volume_nodes not unisolvent with test functions") + raise ValueError("nodes not unisolvent with test functions") test_functions = ( test_basis.derivatives(test_derivative_ax) @@ -439,6 +439,9 @@ def nodal_quadrature_bilinear_form_matrix( if output_nodes is None: output_nodes = input_nodes + if len(test_basis.functions) != output_nodes.shape[1]: + raise ValueError("output_nodes not unisolvent with test functions") + test_functions = ( test_basis.derivatives(test_derivative_ax) if test_derivative_ax is not None else test_basis.functions From d026f656c341a872a7ecbba7b307637d9af7026a Mon Sep 17 00:00:00 2001 From: "Addison A." <70176208+a-alveyblanc@users.noreply.github.com> Date: Wed, 29 Jan 2025 19:27:31 -0600 Subject: [PATCH 24/28] Update modepy/matrices.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- modepy/matrices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index e9a7d153..a6e6c482 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -491,7 +491,7 @@ def spectral_diag_nodal_mass_matrix( return quadrature.weights -# {{{ deprecated remove in 2025-ish +# {{{ deprecated remove in 2026-ish def modal_quad_mass_matrix( quadrature: Quadrature, From b82eea956c9263d7dcf6da9c992dbe85ce9aee86 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Wed, 29 Jan 2025 20:12:49 -0600 Subject: [PATCH 25/28] update according to feedback --- modepy/__init__.py | 4 +- modepy/matrices.py | 80 +++++++++++++++++------------------- modepy/modes.py | 6 ++- modepy/test/test_matrices.py | 10 ++--- modepy/test/test_tools.py | 10 ++--- 5 files changed, 53 insertions(+), 57 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index febbe22c..44ff727a 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -37,7 +37,7 @@ nodal_mass_matrix_for_face, nodal_quad_mass_matrix_for_face, nodal_quadrature_bilinear_form_matrix, - nodal_quadrature_matrix, + nodal_quadrature_test_matrix, resampling_matrix, spectral_diag_nodal_mass_matrix, vandermonde, @@ -164,7 +164,7 @@ "nodal_mass_matrix_for_face", "nodal_quad_mass_matrix_for_face", "nodal_quadrature_bilinear_form_matrix", - "nodal_quadrature_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 e9a7d153..12c64482 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 @@ -76,11 +77,11 @@ """ -BasisFunctionType = Callable[[np.ndarray], np.ndarray] +NodalFunctionType: TypeAlias = Callable[[np.ndarray], np.ndarray] def vandermonde( - functions: Sequence[BasisFunctionType], + functions: Sequence[NodalFunctionType], nodes: np.ndarray ) -> np.ndarray: """Return a (generalized) Vandermonde matrix. @@ -164,7 +165,7 @@ def multi_vandermonde( def resampling_matrix( - basis: Sequence[BasisFunctionType], + basis: Sequence[NodalFunctionType], new_nodes: np.ndarray, old_nodes: np.ndarray, least_squares_ok: bool = False @@ -219,7 +220,7 @@ def is_square(m): def differentiation_matrices( - basis: Sequence[BasisFunctionType], + basis: Sequence[NodalFunctionType], grad_basis: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]], nodes: np.ndarray, from_nodes: np.ndarray | None = None @@ -299,7 +300,7 @@ def diff_matrix_permutation( def inverse_mass_matrix( - basis: Basis | Sequence[BasisFunctionType], + basis: Basis | Sequence[NodalFunctionType], nodes: np.ndarray ) -> np.ndarray: """Return a matrix :math:`A=M^{-1}`, which is the inverse of the one returned @@ -318,7 +319,7 @@ def inverse_mass_matrix( "inverse mass matrix of non-orthogonal basis" ) from None - basis_functions: Sequence[BasisFunctionType] = ( + basis_functions: Sequence[NodalFunctionType] = ( basis.functions) else: basis_functions = basis @@ -334,7 +335,7 @@ def inverse_mass_matrix( def mass_matrix( - basis: Basis | Sequence[BasisFunctionType], + basis: Basis | Sequence[NodalFunctionType], nodes: np.ndarray ) -> np.ndarray: r"""Return a mass matrix :math:`M`, which obeys @@ -352,12 +353,12 @@ def mass_matrix( return la.inv(inverse_mass_matrix(basis, nodes)) -def nodal_quadrature_matrix( +def nodal_quadrature_test_matrix( quadrature: Quadrature, - test_basis: Basis, - test_derivative_ax: int | None = None, + test_basis_functions: Sequence[NodalFunctionType], + test_derivatives: Sequence[NodalFunctionType] | None = None, nodes: np.ndarray | None = None, - mapping_function: Callable[[np.ndarray], np.ndarray] | None = None + mapping_function: NodalFunctionType | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -373,15 +374,15 @@ def nodal_quadrature_matrix( if nodes is None: nodes = quadrature.nodes - if len(test_basis.functions) != nodes.shape[1]: + if len(test_basis_functions) != nodes.shape[1]: raise ValueError("nodes not unisolvent with test functions") test_functions = ( - test_basis.derivatives(test_derivative_ax) - if test_derivative_ax is not None else test_basis.functions + test_derivatives + if test_derivatives is not None else test_basis_functions ) - vdm = vandermonde(test_basis.functions, nodes) + vdm = vandermonde(test_basis_functions, nodes) mapped_nodes = ( mapping_function(quadrature.nodes) @@ -398,13 +399,13 @@ def nodal_quadrature_matrix( def nodal_quadrature_bilinear_form_matrix( quadrature: Quadrature, - test_basis: Basis, - trial_basis: Basis, + test_basis_functions: Sequence[NodalFunctionType], + trial_basis_functions: Sequence[NodalFunctionType], input_nodes: np.ndarray, output_nodes: np.ndarray | None = None, - mapping_function: Callable[[np.ndarray], np.ndarray] | None = None, - test_derivative_ax: int | None = None, - trial_derivative_ax: int | None = None + mapping_function: NodalFunctionType | None = None, + test_derivatives: Sequence[NodalFunctionType] | None = None, + trial_derivatives: Sequence[NodalFunctionType] | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` defined as: @@ -426,30 +427,23 @@ def nodal_quadrature_bilinear_form_matrix( An optional *mapping_function* can be supplied to map the *quadrature* nodes evaluated by *test_basis* before constructing the matrix. - An optional set of *output_nodes* can be supplied if the result should live - on a domain different from the domain of *input_nodes*. - - Test and trial function derivatives can be requested by supplying either - *test_derivative_ax* or *trial_derivative_ax* as an integer corresponding to - an enumerated spatial derivative ax. For example, *test_derivative_ax = 0* - would use *test_basis* differentiated with respect to :math:`x` while - *trial_derivative_ax = 2* would use *trial_basis* differentiated with - respect to :math:`z`. + *test_derivatives* and *trial_derivatives* can be optionally supplied as the + test/trial functions in the inner product. """ if output_nodes is None: output_nodes = input_nodes - if len(test_basis.functions) != output_nodes.shape[1]: + if len(test_basis_functions) != output_nodes.shape[1]: raise ValueError("output_nodes not unisolvent with test functions") test_functions = ( - test_basis.derivatives(test_derivative_ax) - if test_derivative_ax is not None else test_basis.functions + test_derivatives + if test_derivatives is not None else test_basis_functions ) trial_functions = ( - trial_basis.derivatives(trial_derivative_ax) - if trial_derivative_ax is not None else trial_basis.functions + trial_derivatives + if trial_derivatives is not None else trial_basis_functions ) mapped_nodes = ( @@ -464,8 +458,8 @@ def nodal_quadrature_bilinear_form_matrix( quadrature.weights ) - input_vdm = vandermonde(trial_basis.functions, input_nodes) - output_vdm = vandermonde(test_basis.functions, output_nodes) + input_vdm = vandermonde(trial_basis_functions, input_nodes) + output_vdm = vandermonde(test_basis_functions, output_nodes) return la.solve(output_vdm.T, modal_operator @ la.inv(input_vdm)) @@ -491,7 +485,7 @@ def spectral_diag_nodal_mass_matrix( return quadrature.weights -# {{{ deprecated remove in 2025-ish +# {{{ deprecated remove in 2026-ish def modal_quad_mass_matrix( quadrature: Quadrature, @@ -530,8 +524,8 @@ def nodal_quad_mass_matrix( def modal_mass_matrix_for_face( face: Face, face_quad: Quadrature, - trial_functions: Sequence[BasisFunctionType], - test_functions: Sequence[BasisFunctionType] + trial_functions: Sequence[NodalFunctionType], + test_functions: Sequence[NodalFunctionType] ) -> np.ndarray: from warnings import warn warn("`modal_mass_matrix_for_face` is deprecated and will stop working in " @@ -550,8 +544,8 @@ def modal_mass_matrix_for_face( def nodal_mass_matrix_for_face( face: Face, face_quad: Quadrature, - trial_functions: Sequence[BasisFunctionType], - test_functions: Sequence[BasisFunctionType], + trial_functions: Sequence[NodalFunctionType], + test_functions: Sequence[NodalFunctionType], volume_nodes: np.ndarray, face_nodes: np.ndarray ) -> np.ndarray: @@ -576,7 +570,7 @@ def nodal_mass_matrix_for_face( def modal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[BasisFunctionType], + test_functions: Sequence[NodalFunctionType], ) -> np.ndarray: from warnings import warn warn("`modal_quad_mass_matrix_for_face` is deprecated and will stop working " @@ -593,7 +587,7 @@ def modal_quad_mass_matrix_for_face( def nodal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[BasisFunctionType], + test_functions: Sequence[NodalFunctionType], volume_nodes: np.ndarray, ) -> np.ndarray: from warnings import warn diff --git a/modepy/modes.py b/modepy/modes.py index 05a31349..917da1c5 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 ca130925..40b3ef90 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -114,10 +114,10 @@ def test_bilinear_forms( weak_operator = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, - test_basis=basis, - trial_basis=basis, + test_basis_functions=basis.functions, + trial_basis_functions=basis.functions, input_nodes=nodes, - test_derivative_ax=ax + test_derivatives=basis.derivatives(ax) ) err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) @@ -125,8 +125,8 @@ def test_bilinear_forms( else: quad_mass_mat = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, - test_basis=basis, - trial_basis=basis, + test_basis_functions=basis.functions, + trial_basis_functions=basis.functions, input_nodes=nodes ) diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py index 8dc8d437..b298ece0 100644 --- a/modepy/test/test_tools.py +++ b/modepy/test/test_tools.py @@ -304,7 +304,7 @@ def test_nodal_quadrature_bilinear_form_matrix_for_face(dims, shape_cls, order=3 from modepy.matrices import ( nodal_quadrature_bilinear_form_matrix, - nodal_quadrature_matrix, + nodal_quadrature_test_matrix, ) for face in mp.faces_for_shape(vol_shape): face_space = mp.space_for_shape(face, order) @@ -316,16 +316,16 @@ def test_nodal_quadrature_bilinear_form_matrix_for_face(dims, shape_cls, order=3 fmm = nodal_quadrature_bilinear_form_matrix( quadrature=face_quad, - test_basis=volume_basis, - trial_basis=face_basis, + test_basis_functions=volume_basis.functions, + trial_basis_functions=face_basis.functions, input_nodes=face_nodes, output_nodes=volume_nodes, mapping_function=face.map_to_volume ) - fmm2 = nodal_quadrature_matrix( + fmm2 = nodal_quadrature_test_matrix( quadrature=face_quad2, - test_basis=volume_basis, + test_basis_functions=volume_basis.functions, nodes=volume_nodes, mapping_function=face.map_to_volume ) From 461708813ac15635319eba7d35610d2e333462bc Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Wed, 29 Jan 2025 20:17:22 -0600 Subject: [PATCH 26/28] add back quad_mass_matrix test --- modepy/__init__.py | 2 ++ modepy/test/test_matrices.py | 43 ++++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/modepy/__init__.py b/modepy/__init__.py index 44ff727a..c4594899 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -35,6 +35,7 @@ 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, @@ -162,6 +163,7 @@ "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", diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index 40b3ef90..ce7fb3a1 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -33,6 +33,49 @@ import modepy as mp +@pytest.mark.parametrize("shape", [ + mp.Simplex(1), + mp.Simplex(2), + mp.Simplex(3), + mp.Hypercube(2), + mp.Hypercube(3), + ]) +@pytest.mark.parametrize("order", [0, 1, 2, 4]) +@pytest.mark.parametrize("nodes_on_bdry", [False, True]) +def test_nodal_mass_matrix_against_quad( + shape: mp.Shape, order: int, nodes_on_bdry: bool, + ) -> None: + 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 isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1): + if nodes_on_bdry: + nodes = mp.legendre_gauss_lobatto_tensor_product_nodes( + shape.dim, order, + ) + else: + nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order) + elif isinstance(shape, mp.Simplex): + if nodes_on_bdry: + nodes = mp.warp_and_blend_nodes(shape.dim, order) + else: + nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes + + else: + raise AssertionError() + + basis = mp.orthonormal_basis_for_space(space, shape) + + quad_mass_mat = mp.nodal_quad_mass_matrix(quad, basis.functions, nodes) + vdm_mass_mat = mp.mass_matrix(basis, nodes) + + nodes_to_quad = mp.resampling_matrix(basis.functions, quad.nodes, nodes) + + err = la.norm(quad_mass_mat@nodes_to_quad - vdm_mass_mat)/la.norm(vdm_mass_mat) + assert err < 1e-14 + + @pytest.mark.parametrize("shape", [ mp.Hypercube(1), mp.Hypercube(2), From d6d53fad0f262a4e42a6efcbfce07f9a1cdfa5f8 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Wed, 29 Jan 2025 20:21:43 -0600 Subject: [PATCH 27/28] fix docs --- modepy/matrices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 12c64482..2613c404 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -63,7 +63,7 @@ .. autofunction:: inverse_mass_matrix .. autofunction:: mass_matrix .. autofunction:: spectral_diag_nodal_mass_matrix -.. autofunction:: nodal_quadrature_matrix +.. autofunction:: nodal_quadrature_test_matrix .. autofunction:: nodal_quadrature_bilinear_form_matrix Differentiation is also convenient to express by using :math:`V^{-1}` to From 19206fc1525179213734f0cf849e446ac5e384fa Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Tue, 4 Feb 2025 23:54:19 -0600 Subject: [PATCH 28/28] improve arg names and docs --- modepy/matrices.py | 105 ++++++++++++++++------------------- modepy/test/test_matrices.py | 13 +++-- modepy/test/test_tools.py | 13 +++-- 3 files changed, 64 insertions(+), 67 deletions(-) diff --git a/modepy/matrices.py b/modepy/matrices.py index 2613c404..2d3da118 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -77,11 +77,11 @@ """ -NodalFunctionType: TypeAlias = Callable[[np.ndarray], np.ndarray] +NodalFunction: TypeAlias = Callable[[np.ndarray], np.ndarray] def vandermonde( - functions: Sequence[NodalFunctionType], + functions: Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: """Return a (generalized) Vandermonde matrix. @@ -165,7 +165,7 @@ def multi_vandermonde( def resampling_matrix( - basis: Sequence[NodalFunctionType], + basis: Sequence[NodalFunction], new_nodes: np.ndarray, old_nodes: np.ndarray, least_squares_ok: bool = False @@ -220,7 +220,7 @@ def is_square(m): def differentiation_matrices( - basis: Sequence[NodalFunctionType], + basis: Sequence[NodalFunction], grad_basis: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]], nodes: np.ndarray, from_nodes: np.ndarray | None = None @@ -300,7 +300,7 @@ def diff_matrix_permutation( def inverse_mass_matrix( - basis: Basis | Sequence[NodalFunctionType], + basis: Basis | Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: """Return a matrix :math:`A=M^{-1}`, which is the inverse of the one returned @@ -319,7 +319,7 @@ def inverse_mass_matrix( "inverse mass matrix of non-orthogonal basis" ) from None - basis_functions: Sequence[NodalFunctionType] = ( + basis_functions: Sequence[NodalFunction] = ( basis.functions) else: basis_functions = basis @@ -335,7 +335,7 @@ def inverse_mass_matrix( def mass_matrix( - basis: Basis | Sequence[NodalFunctionType], + basis: Basis | Sequence[NodalFunction], nodes: np.ndarray ) -> np.ndarray: r"""Return a mass matrix :math:`M`, which obeys @@ -355,10 +355,10 @@ def mass_matrix( def nodal_quadrature_test_matrix( quadrature: Quadrature, - test_basis_functions: Sequence[NodalFunctionType], - test_derivatives: Sequence[NodalFunctionType] | None = None, + test_functions: Sequence[NodalFunction], + nodal_interp_functions: Sequence[NodalFunction], nodes: np.ndarray | None = None, - mapping_function: NodalFunctionType | None = None + test_function_node_map: NodalFunction | None = None ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` that satisfies: @@ -370,42 +370,39 @@ def nodal_quadrature_test_matrix( *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. + + *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 len(test_basis_functions) != nodes.shape[1]: - raise ValueError("nodes not unisolvent with test functions") - - test_functions = ( - test_derivatives - if test_derivatives is not None else test_basis_functions - ) + if len(nodal_interp_functions) != nodes.shape[1]: + raise ValueError("nodes not unisolvent with nodal_interp_functions") - vdm = vandermonde(test_basis_functions, nodes) + vdm = vandermonde(test_functions, nodes) - mapped_nodes = ( - mapping_function(quadrature.nodes) - if mapping_function is not None else quadrature.nodes + test_nodes = ( + test_function_node_map(quadrature.nodes) + if test_function_node_map is not None else quadrature.nodes ) - modal_operator = np.array([ - test_f(mapped_nodes) * quadrature.weights - for test_f in test_functions - ]) + modal_mat = vandermonde(test_functions, test_nodes).T*quadrature.weights - return la.solve(vdm.T, modal_operator) + return la.solve(vdm.T, modal_mat) def nodal_quadrature_bilinear_form_matrix( quadrature: Quadrature, - test_basis_functions: Sequence[NodalFunctionType], - trial_basis_functions: Sequence[NodalFunctionType], + 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, - mapping_function: NodalFunctionType | None = None, - test_derivatives: Sequence[NodalFunctionType] | None = None, - trial_derivatives: Sequence[NodalFunctionType] | None = None + test_function_node_map: NodalFunction | None = None, ) -> np.ndarray: r"""Using *quadrature*, provide a matrix :math:`A` defined as: @@ -424,31 +421,25 @@ def nodal_quadrature_bilinear_form_matrix( where :math:`u_i` are nodal coefficients. - An optional *mapping_function* can be supplied to map the *quadrature* nodes - evaluated by *test_basis* before constructing the matrix. - - *test_derivatives* and *trial_derivatives* can be optionally supplied as the - test/trial functions in the inner product. + *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 output_nodes is None: output_nodes = input_nodes - if len(test_basis_functions) != output_nodes.shape[1]: - raise ValueError("output_nodes not unisolvent with test functions") - - test_functions = ( - test_derivatives - if test_derivatives is not None else test_basis_functions - ) + if len(nodal_interp_functions_test) != output_nodes.shape[1]: + raise ValueError( + "output_nodes not unisolvent with nodal_test_interp_functions") - trial_functions = ( - trial_derivatives - if trial_derivatives is not None else trial_basis_functions - ) + if len(nodal_interp_functions_trial) != input_nodes.shape[1]: + raise ValueError( + "input_nodes not unisolvent with nodal_trial_interp_functions") mapped_nodes = ( - mapping_function(quadrature.nodes) - if mapping_function is not None else quadrature.nodes + test_function_node_map(quadrature.nodes) + if test_function_node_map is not None else quadrature.nodes ) modal_operator = np.einsum( @@ -458,8 +449,8 @@ def nodal_quadrature_bilinear_form_matrix( quadrature.weights ) - input_vdm = vandermonde(trial_basis_functions, input_nodes) - output_vdm = vandermonde(test_basis_functions, output_nodes) + 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)) @@ -524,8 +515,8 @@ def nodal_quad_mass_matrix( def modal_mass_matrix_for_face( face: Face, face_quad: Quadrature, - trial_functions: Sequence[NodalFunctionType], - test_functions: Sequence[NodalFunctionType] + 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 " @@ -544,8 +535,8 @@ def modal_mass_matrix_for_face( def nodal_mass_matrix_for_face( face: Face, face_quad: Quadrature, - trial_functions: Sequence[NodalFunctionType], - test_functions: Sequence[NodalFunctionType], + trial_functions: Sequence[NodalFunction], + test_functions: Sequence[NodalFunction], volume_nodes: np.ndarray, face_nodes: np.ndarray ) -> np.ndarray: @@ -570,7 +561,7 @@ def nodal_mass_matrix_for_face( def modal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[NodalFunctionType], + test_functions: Sequence[NodalFunction], ) -> np.ndarray: from warnings import warn warn("`modal_quad_mass_matrix_for_face` is deprecated and will stop working " @@ -587,7 +578,7 @@ def modal_quad_mass_matrix_for_face( def nodal_quad_mass_matrix_for_face( face: Face, face_quad: Quadrature, - test_functions: Sequence[NodalFunctionType], + test_functions: Sequence[NodalFunction], volume_nodes: np.ndarray, ) -> np.ndarray: from warnings import warn diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py index ce7fb3a1..6e326fe1 100644 --- a/modepy/test/test_matrices.py +++ b/modepy/test/test_matrices.py @@ -157,10 +157,11 @@ def test_bilinear_forms( weak_operator = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, - test_basis_functions=basis.functions, - trial_basis_functions=basis.functions, + test_functions=basis.derivatives(ax), + trial_functions=basis.functions, + nodal_interp_functions_test=basis.functions, + nodal_interp_functions_trial=basis.functions, input_nodes=nodes, - test_derivatives=basis.derivatives(ax) ) err = la.norm(mass_inv @ weak_operator.T @ f - fp) / la.norm(fp) @@ -168,8 +169,10 @@ def test_bilinear_forms( else: quad_mass_mat = mp.nodal_quadrature_bilinear_form_matrix( quadrature=quad, - test_basis_functions=basis.functions, - trial_basis_functions=basis.functions, + test_functions=basis.functions, + trial_functions=basis.functions, + nodal_interp_functions_test=basis.functions, + nodal_interp_functions_trial=basis.functions, input_nodes=nodes ) diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py index b298ece0..2d22f8a9 100644 --- a/modepy/test/test_tools.py +++ b/modepy/test/test_tools.py @@ -316,18 +316,21 @@ def test_nodal_quadrature_bilinear_form_matrix_for_face(dims, shape_cls, order=3 fmm = nodal_quadrature_bilinear_form_matrix( quadrature=face_quad, - test_basis_functions=volume_basis.functions, - trial_basis_functions=face_basis.functions, + 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, - mapping_function=face.map_to_volume + test_function_node_map=face.map_to_volume ) fmm2 = nodal_quadrature_test_matrix( quadrature=face_quad2, - test_basis_functions=volume_basis.functions, + test_functions=volume_basis.functions, + nodal_interp_functions=volume_basis.functions, nodes=volume_nodes, - mapping_function=face.map_to_volume + test_function_node_map=face.map_to_volume ) for f_face in face_basis.functions: