diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index 99557843..987cf416 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -435,6 +435,26 @@ template T FunctionTree::integrate() const { return jacobian * result; } +template T FunctionTree::integrate(int dim, bool largerSide) const { + T result = 0.0; + for (int i = 0; i < this->rootBox.size(); i++) { + const FunctionNode &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 diff --git a/src/trees/FunctionTree.h b/src/trees/FunctionTree.h index 9d976d6b..ffc7cd45 100644 --- a/src/trees/FunctionTree.h +++ b/src/trees/FunctionTree.h @@ -62,6 +62,15 @@ template class FunctionTree final : public MWTree, 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 &r); T evalf(const Coord &r) const override; diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 250796c2..ab374608 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -187,6 +187,15 @@ template ComplexDouble CompFunction::integrate() const { return integral; } +template ComplexDouble CompFunction::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 double CompFunction::norm() const { double norm = getSquareNorm(); if (norm > 0.0) norm = std::sqrt(norm); diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 3b63a63b..3979866a 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -126,6 +126,15 @@ template 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); diff --git a/tests/trees/function_tree.cpp b/tests/trees/function_tree.cpp index 99e5e709..ab657f8f 100644 --- a/tests/trees/function_tree.cpp +++ b/tests/trees/function_tree.cpp @@ -57,6 +57,7 @@ template 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);