Add generic nodal/modal bilinear form routines

modepy/__init__.py

@@ -32,12 +32,16 @@
     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_quadrature_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,
     resampling_matrix,
     spectral_diag_nodal_mass_matrix,
     vandermonde,
@@ -158,13 +162,17 @@
     "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_quadrature_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",
     "node_tuples_for_space",
     "orthonormal_basis_for_space",
     "quadrature_for_space",

modepy/matrices.py

@@ -66,6 +66,10 @@
 .. 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_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
@@ -351,6 +355,117 @@ def mass_matrix(
     return la.inv(inverse_mass_matrix(basis, nodes))
 
 
+def modal_quadrature_operator(
+        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`, 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_quadrature_operator(
+        quadrature: Quadrature,
+        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:
+
+    .. 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
+    solution at the *quadrature* nodes.
+    """
+    if nodes is None:
+        nodes = quadrature.nodes
+
+    if len(proj_functions) != nodes.shape[1]:
+        raise ValueError("volume_nodes not unisolvent with proj_functions")
+
+    vdm = vandermonde(proj_functions, nodes)
+
+    return la.solve(
+        vdm.T, modal_quadrature_operator(
+            quadrature,
+            test_functions
+        )
+    )
+
+
+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 \psi_i(r_k)
+        \phi_j(r_k)) u_j,
+
+    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 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*.
+    """
+    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]],
@@ -367,12 +482,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(
@@ -396,16 +506,8 @@ def nodal_quad_mass_matrix(
 
     .. versionadded :: 2024.2
     """
-    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_quadrature_operator(
+        quadrature, test_functions, test_functions, nodes)
 
 
 def spectral_diag_nodal_mass_matrix(
@@ -549,5 +651,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))
 
-
 # vim: foldmethod=marker

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)
 
@@ -118,6 +112,65 @@ def test_tensor_product_diag_mass_matrix(shape: mp.Shape) -> None:
                 assert np.max(err) < 1e-14
 
 
+@pytest.mark.parametrize("shape_cls", [mp.Hypercube, mp.Simplex])
+@pytest.mark.parametrize("dim", [1, 2, 3])
+@pytest.mark.parametrize("order", [0, 1, 2, 4])
+@pytest.mark.parametrize("nodes_on_bdry", [False, True])
+@pytest.mark.parametrize("test_weak_d_dr", [False, True])
+def test_bilinear_forms(
+            shape_cls: type[mp.Shape],
+            dim: int,
+            order: int,
+            nodes_on_bdry: bool,
+            test_weak_d_dr: bool
+        ) -> None:
+    shape = shape_cls(dim)
+    space = mp.space_for_shape(shape, order)
+
+    quad_space = mp.space_for_shape(shape, 2*order)
+    quad = mp.quadrature_for_space(quad_space, shape)
+
+    if nodes_on_bdry:
+        nodes = mp.edge_clustered_nodes_for_space(space, shape)
+    else:
+        if isinstance(shape, mp.Hypercube) or shape == mp.Simplex(1):
+            nodes = mp.legendre_gauss_tensor_product_nodes(shape.dim, order)
+        elif isinstance(shape, mp.Simplex):
+            nodes = mp.VioreanuRokhlinSimplexQuadrature(order, shape.dim).nodes
+        else:
+            raise AssertionError()
+
+    basis = mp.orthonormal_basis_for_space(space, shape)
+
+    if test_weak_d_dr and order not in [0, 1]:
+        mass_inv = mp.inverse_mass_matrix(basis, nodes)
+
+        for ax in range(dim):
+            f = 1 - nodes[ax]**2
+            fp = -2*nodes[ax]
+
+            weak_operator = mp.nodal_quad_bilinear_form(
+                quad,
+                basis.derivatives(ax),
+                basis.functions,
+                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(
+            quad,
+            basis.functions,
+            basis.functions,
+            basis.functions,
+            nodes)
+
+        vdm_mass_mat = mp.mass_matrix(basis, nodes)
+        err = la.norm(quad_mass_mat - vdm_mass_mat) / la.norm(vdm_mass_mat)
+        assert err < 1e-14
+
+
 # You can test individual routines by typing
 # $ python test_modes.py 'test_routine()'