Expose reference moment fitting systems to Python

Author: Q-Minh

bindings/pypbat/math/MomentFitting.cpp

@@ -52,6 +52,53 @@ void BindMomentFitting(pybind11::module& m)
         "of X1. S2 is the same for columns of X2. w2 are the quadrature weights "
         "associated with X2. If with_error==True, the polynomial integration error "
         "is computed and returned.");
+
+    m.def(
+        "reference_moment_fitting_systems",
+        [](Eigen::Ref<IndexVectorX const> const& S1,
+           Eigen::Ref<MatrixX const> const& X1,
+           Eigen::Ref<IndexVectorX const> const& S2,
+           Eigen::Ref<MatrixX const> const& X2,
+           Eigen::Ref<VectorX const> const& w2,
+           Index order) {
+            MatrixX M, B;
+            IndexVectorX P;
+            common::ForRange<1, 4>([&]<auto Order>() {
+                if (order == Order)
+                {
+                    std::tie(M, B, P) =
+                        pbat::math::ReferenceMomentFittingSystems<Order>(S1, X1, S2, X2, w2);
+                }
+            });
+            if (M.size() == 0)
+                throw std::invalid_argument("transfer_quadrature only accepts 1 <= order <= 4.");
+            return std::make_tuple(M, B, P);
+        },
+        pyb::arg("S1"),
+        pyb::arg("X1"),
+        pyb::arg("S2"),
+        pyb::arg("X2"),
+        pyb::arg("w2"),
+        pyb::arg("order") = 1,
+        "Obtain a collection of reference moment fitting systems (M, B, P), where M[:, "
+        "P[s]:P[s+1]] is the reference moment fitting matrix for simplex s, and b[:,s] is its "
+        "corresponding right-hand side. X1, S1 are the |#dims|x|#quad.pts.| array of quadrature "
+        "points and corresponding simplices, and X2,w2,S2 are the |#dims|x|#old quad.pts.| array "
+        "of existing quadrature points and corresponding simplices, with weights w2. order "
+        "specifies the polynomial order of integration that we wish to reproduce exactly.");
+
+    m.def(
+        "block_diagonalize_moment_fitting",
+        [](Eigen::Ref<MatrixX const> const& M,
+           Eigen::Ref<MatrixX const> const& B,
+           Eigen::Ref<IndexVectorX const> const& P) {
+            return pbat::math::BlockDiagonalReferenceMomentFittingSystem(M, B, P);
+        },
+        pyb::arg("M"),
+        pyb::arg("B"),
+        pyb::arg("P"),
+        "Assemble the block diagonal row sparse matrix GM, such that GM @ w = B.flatten(order='F') "
+        "contains all the reference moment fitting systems in (M,B,P).");
 }
 
 } // namespace math