Population analysis

src/trees/FunctionTree.cpp

@@ -435,6 +435,26 @@ template <int D, typename T> T FunctionTree<D, T>::integrate() const {
     return jacobian * result;
 }
 
+template <int D, typename T> T FunctionTree<D, T>::integrate(int dim, bool largerSide) const {
+    T result = 0.0;
+    for (int i = 0; i < this->rootBox.size(); i++) {
+        const FunctionNode<D, T> &fNode = getRootFuncNode(i);
+        if (fNode.getLowerBounds()[dim] == this->rootBox.getLowerBound(dim) && not largerSide)
+            result += fNode.integrate();
+        if (fNode.getUpperBounds()[dim] == this->rootBox.getUpperBound(dim) && largerSide)
+            result += fNode.integrate();
+    }
+
+    // Handle potential scaling
+    auto scaling_factor = this->getMRA().getWorldBox().getScalingFactors();
+    auto jacobian = 1.0;
+    for (const auto &sf_i : scaling_factor) jacobian *= std::sqrt(sf_i);
+    // Square root of scaling factor in each diection. The seemingly missing
+    // multiplication by the square root of sf_i is included in the basis
+
+    return jacobian * result;
+}
+
 /** @returns Integral of a representable function over the grid given by the tree */
 template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunction_M &f) {
     // traverse tree, and treat end nodes only

src/trees/FunctionTree.h

@@ -62,6 +62,15 @@ template <int D, typename T> class FunctionTree final : public MWTree<D, T>, pub
     ~FunctionTree() override;
 
     T integrate() const;
+
+    /** 
+     * @brief Integrate over half of the space
+     * @param dim Dimension along which to split
+     * @param largerSide If true, integrate over the larger side (x>0 if dim=0)
+     * @return Integral of the function over parts of the computational domain
+     */
+    T integrate(int dim, bool largerSide) const;
+
     double integrateEndNodes(RepresentableFunction_M &f);
     T evalf_precise(const Coord<D> &r);
     T evalf(const Coord<D> &r) const override;

src/utils/CompFunction.cpp

@@ -187,6 +187,15 @@ template <int D> ComplexDouble CompFunction<D>::integrate() const {
     return integral;
 }
 
+template <int D> ComplexDouble CompFunction<D>::integrate(int dim, bool largerSide) const {
+    ComplexDouble integral;
+    if (isreal())
+        integral = CompD[0]->integrate(dim, largerSide);
+    else
+        integral = CompC[0]->integrate(dim, largerSide);
+    return integral;
+}
+
 template <int D> double CompFunction<D>::norm() const {
     double norm = getSquareNorm();
     if (norm > 0.0) norm = std::sqrt(norm);

src/utils/CompFunction.h

@@ -126,6 +126,15 @@ template <int D> class CompFunction {
 
     CompFunction paramCopy(bool alloc = false) const;
     ComplexDouble integrate() const;
+
+    /**
+     * @brief Integrate over half of the space
+     * @param dim Dimension along which to split
+     * @param largerSide If true, integrate over the larger side (x>0 if dim=0)
+     * @return Resulting integral
+     */
+    ComplexDouble integrate(int dim, bool largerSide) const;
+
     double norm() const;
     double getSquareNorm() const;
     void alloc(int nalloc = 1, bool zero = true);

tests/trees/function_tree.cpp

@@ -57,6 +57,7 @@ template <int D> void testZeroFunction() {
         }
         THEN("its squared norm is zero") { REQUIRE(tree.getSquareNorm() == Catch::Approx(0.0)); }
         THEN("it integrates to zero") { REQUIRE(tree.integrate() == Catch::Approx(0.0)); }
+        THEN("half of it integrates to zero") { REQUIRE(tree.integrate(0, true) == Catch::Approx(0.0)); }
         THEN("the dot product with itself is zero") { REQUIRE(dot(tree, tree) == Catch::Approx(0.0)); }
     }
     finalize(&mra);