Expose quadrature weight transfer to Python

Author: Q-Minh

bindings/pypbat/math/CMakeLists.txt

@@ -3,11 +3,13 @@ target_sources(PhysicsBasedAnimationToolkit_Python
     FILE_SET api
     FILES
     "Math.h"
+    "MomentFitting.h"
 )
 
 target_sources(PhysicsBasedAnimationToolkit_Python
     PRIVATE
     "Math.cpp"
+    "MomentFitting.cpp"
 )
 
 add_subdirectory(linalg)

bindings/pypbat/math/Math.cpp

@@ -1,5 +1,6 @@
 #include "Math.h"
 
+#include "MomentFitting.h"
 #include "linalg/LinAlg.h"
 
 #include <string>
@@ -10,7 +11,7 @@ namespace math {
 
 void Bind(pybind11::module& m)
 {
-    namespace pyb = pybind11;
+    BindMomentFitting(m);
     auto mlinalg = m.def_submodule("linalg");
     linalg::Bind(mlinalg);
 }

bindings/pypbat/math/MomentFitting.cpp

@@ -0,0 +1,59 @@
+#include "MomentFitting.h"
+
+#include <exception>
+#include <pbat/common/ConstexprFor.h>
+#include <pbat/math/MomentFitting.h>
+#include <pybind11/eigen.h>
+
+namespace pbat {
+namespace py {
+namespace math {
+
+void BindMomentFitting(pybind11::module& m)
+{
+    namespace pyb = pybind11;
+    m.def(
+        "transfer_quadrature",
+        [](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,
+           bool bEvaluateError,
+           Index maxIterations,
+           Scalar precision) {
+            VectorX w1, err;
+            common::ForRange<1, 4>([&]<auto Order>() {
+                if (order == Order)
+                {
+                    std::tie(w1, err) = pbat::math::TransferQuadrature<
+                        Order>(S1, X1, S2, X2, w2, bEvaluateError, maxIterations, precision);
+                }
+            });
+            if (w1.size() == 0)
+                throw std::invalid_argument("transfer_quadrature only accepts 1 <= order <= 4.");
+            return std::make_pair(w1, err);
+        },
+        pyb::arg("S1"),
+        pyb::arg("X1"),
+        pyb::arg("S2"),
+        pyb::arg("X2"),
+        pyb::arg("w2"),
+        pyb::arg("order")      = 1,
+        pyb::arg("with_error") = false,
+        pyb::arg("max_iters")  = 20,
+        pyb::arg("precision")  = std::numeric_limits<Scalar>::epsilon(),
+        "Obtain weights w1 by transferring an existing quadrature rule (X2,w2) "
+        "defined on a domain composed of simplices onto a new quadrature rule "
+        "(X1,w1) defined on the same domain, given fixed quadrature points X1. "
+        "S1 is a |X1.shape[1]| array of simplex indices associated with the "
+        "corresponding quadrature point in columns (i.e. the quadrature points) "
+        "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.");
+}
+
+} // namespace math
+} // namespace py
+} // namespace pbat

bindings/pypbat/math/MomentFitting.h

@@ -0,0 +1,16 @@
+#ifndef PYPBAT_MATH_MOMENT_FITTING_H
+#define PYPBAT_MATH_MOMENT_FITTING_H
+
+#include <pybind11/pybind11.h>
+
+namespace pbat {
+namespace py {
+namespace math {
+
+void BindMomentFitting(pybind11::module& m);
+
+} // namespace math
+} // namespace py
+} // namespace pbat
+
+#endif // PYPBAT_MATH_MOMENT_FITTING_H

python/pbatoolkit/math/__init__.py

@@ -0,0 +1,8 @@
+from .._pbat import math as _math
+import sys
+import inspect
+
+__module = sys.modules[__name__]
+for _name, _attr in inspect.getmembers(_math):
+    if not _name.startswith("__"):
+        setattr(__module, _name, _attr)