Fix quadrature tests and add and order evaluation functionality

tests/test_FunctionEvaluation.cpp

@@ -38,16 +38,29 @@ TEST_CASE("QuadratureFormula", "[analysis]")
       auto point_list  = quadrature_formula.pointList();
       auto weight_list = quadrature_formula.weightList();
 
-      auto x = [&a, &b](auto x_hat) { return 0.5 * (b - a) * (x_hat + 1); };
+      double alpha = 0.5 * (b - a);
+      double beta  = 0.5 * (a + b);
+
+      auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
 
       auto value = weight_list[0] * f(x(point_list[0]));
       for (size_t i = 1; i < weight_list.size(); ++i) {
         value += weight_list[i] * f(x(point_list[i]));
       }
-      return 0.5 * (b - a) * value;
+
+      return alpha * value;
+    };
+
+    auto get_order = [&integrate](auto f, auto quadrature_formula, const double a, const double b,
+                                  const double exact_value) {
+      double int_ab          = integrate(f, quadrature_formula, a, b);
+      double int_first_half  = integrate(f, quadrature_formula, a, 0.5 * (a + b));
+      double int_second_half = integrate(f, quadrature_formula, 0.5 * (a + b), b);
+
+      return -std::log((int_first_half + int_second_half - exact_value) / (int_ab - exact_value)) / std::log(2);
     };
 
-    auto p0 = [](const TinyVector<1>&) { return 3; };
+    auto p0 = [](const TinyVector<1>&) { return 1; };
     auto p1 = [](const TinyVector<1>& X) {
       const double x = X[0];
       return 2 * x + 1;
@@ -56,6 +69,14 @@ TEST_CASE("QuadratureFormula", "[analysis]")
       const double x = X[0];
       return 3 * x * x + 2 * x + 1;
     };
+    auto p3 = [](const TinyVector<1>& X) {
+      const double x = X[0];
+      return 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1;
+    };
+    auto p4 = [](const TinyVector<1>& X) {
+      const double x = X[0];
+      return 5 * std::pow(x, 4) + 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1;
+    };
 
     SECTION("exact integration")
     {
@@ -67,8 +88,11 @@ TEST_CASE("QuadratureFormula", "[analysis]")
         REQUIRE(l0.numberOfPoints() == 1);
         REQUIRE(l1.numberOfPoints() == 1);
 
-        REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(1.5));
+        REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(0.5));
         REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
+        REQUIRE(integrate(p2, l1, 0, 1) != Catch::Approx(3));
+
+        REQUIRE(get_order(p2, l1, 0, 1, 3) == Catch::Approx(2));
       }
 
       SECTION("order 2 and 3")
@@ -79,9 +103,13 @@ TEST_CASE("QuadratureFormula", "[analysis]")
         REQUIRE(l2.numberOfPoints() == 2);
         REQUIRE(l3.numberOfPoints() == 2);
 
-        REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(1.5));
+        REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(0.5));
         REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
         REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
+        REQUIRE(integrate(p3, l3, 0, 1) == Catch::Approx(4));
+        REQUIRE(integrate(p4, l3, 0, 1) != Catch::Approx(5));
+
+        REQUIRE(get_order(p4, l3, 0, 1, 5) == Catch::Approx(4));
       }
     }
   }