diff --git a/tests/test_EdgeIntegrator.cpp b/tests/test_EdgeIntegrator.cpp
index f8ea09cee6bdeef7352281f207e2dfcc3b66e408..5dfccbc6c7427a2bd27634c49957209c61c6f3ce 100644
--- a/tests/test_EdgeIntegrator.cpp
+++ b/tests/test_EdgeIntegrator.cpp
@@ -1,6 +1,7 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
+#include <algebra/TinyMatrix.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
@@ -15,242 +16,745 @@
TEST_CASE("EdgeIntegrator", "[scheme]")
{
- SECTION("3D")
+ SECTION("scalar")
{
- using R3 = TinyVector<3>;
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ auto f = [](const R3& X) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
+ };
+
+ std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
+ mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
+ mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+
+ for (auto mesh_info : mesh_list) {
+ auto mesh_name = mesh_info.first;
+ auto mesh = mesh_info.second;
+
+ Array<const double> int_f_per_edge = [=] {
+ Array<double> int_f(mesh->numberOfEdges());
+ auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+ auto edge_node_list = edge_to_node_matrix[edge_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]);
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.velocityNorm() * f(T(xi));
+ }
- auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+ int_f[edge_id] = integral;
+ });
- auto f = [](const R3& X) -> double {
- const double x = X[0];
- const double y = X[1];
- const double z = X[2];
- return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
- };
+ return int_f;
+ }();
- std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
- mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
- mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
-
- for (auto mesh_info : mesh_list) {
- auto mesh_name = mesh_info.first;
- auto mesh = mesh_info.second;
-
- Array<const double> int_f_per_edge = [=] {
- Array<double> int_f(mesh->numberOfEdges());
- auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix();
-
- parallel_for(
- mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
- auto edge_node_list = edge_to_node_matrix[edge_id];
- auto xr = mesh->xr();
- double integral = 0;
-
- LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]);
- auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.velocityNorm() * f(T(xi));
- }
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
+ {
+ SECTION("all edges")
+ {
+ SECTION("EdgeValue")
+ {
+ EdgeValue<double> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
- int_f[edge_id] = integral;
- });
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
- return int_f;
- }();
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- SECTION(mesh_name)
- {
- SECTION("direct formula")
- {
- SECTION("all edges")
- {
- SECTION("EdgeValue")
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfEdges());
+
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("edge list")
{
- EdgeValue<double> values(mesh->connectivity());
- EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
- values);
+ SECTION("Array")
+ {
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ REQUIRE(k == edge_list.size());
+ }
+
+ Array<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<double> values =
+ EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
+ }
- SECTION("Array")
+ SECTION("tensorial formula")
+ {
+ SECTION("all edges")
{
- Array<double> values(mesh->numberOfEdges());
+ SECTION("EdgeValue")
+ {
+ EdgeValue<double> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
- EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfEdges());
+
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("edge list")
{
- SmallArray<double> values(mesh->numberOfEdges());
- EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ SECTION("Array")
+ {
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ REQUIRE(k == edge_list.size());
+ }
+
+ Array<double> values =
+ EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<double> values =
+ EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
+ }
+ }
+ }
+ }
+
+ SECTION("vector")
+ {
+ using R2 = TinyVector<2>;
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ auto f = [](const R3& X) -> R2 {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return R2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x + 3 * y - 1};
+ };
+
+ std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
+ mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
+ mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+
+ for (auto mesh_info : mesh_list) {
+ auto mesh_name = mesh_info.first;
+ auto mesh = mesh_info.second;
+
+ Array<const R2> int_f_per_edge = [=] {
+ Array<R2> int_f(mesh->numberOfEdges());
+ auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+ auto edge_node_list = edge_to_node_matrix[edge_id];
+ auto xr = mesh->xr();
+ R2 integral = zero;
+
+ LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]);
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.velocityNorm() * f(T(xi));
+ }
+
+ int_f[edge_id] = integral;
+ });
- SECTION("edge list")
+ return int_f;
+ }();
+
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
{
- SECTION("Array")
+ SECTION("all edges")
{
- Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+ SECTION("EdgeValue")
+ {
+ EdgeValue<R2> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<R2> values(mesh->numberOfEdges());
+
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ SECTION("SmallArray")
{
- size_t k = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
- edge_list[k] = edge_id;
+ SmallArray<R2> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
}
- REQUIRE(k == edge_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
+ }
- Array<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+ SECTION("edge list")
+ {
+ SECTION("Array")
+ {
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ Array<R2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
- double error = 0;
- for (size_t i = 0; i < edge_list.size(); ++i) {
- error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<R2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
+ }
- SECTION("SmallArray")
+ SECTION("tensorial formula")
+ {
+ SECTION("all edges")
{
- SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+ SECTION("EdgeValue")
+ {
+ EdgeValue<R2> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
{
- size_t k = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
- edge_list[k] = edge_id;
+ Array<R2> values(mesh->numberOfEdges());
+
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
}
- REQUIRE(k == edge_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- SmallArray<double> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+ SECTION("SmallArray")
+ {
+ SmallArray<R2> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < edge_list.size(); ++i) {
- error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
+ }
+
+ SECTION("edge list")
+ {
+ SECTION("Array")
+ {
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ REQUIRE(k == edge_list.size());
+ }
+
+ Array<R2> values = EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<R2> values =
+ EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
}
+ }
+ }
+ }
- SECTION("tensorial formula")
+ SECTION("matrix")
+ {
+ using R2x2 = TinyMatrix<2>;
+ auto l2Norm = [](const R2x2& A) -> double {
+ return std::sqrt(A(0, 0) * A(0, 0) + A(1, 0) * A(1, 0) + A(0, 1) * A(0, 1) + A(1, 1) * A(1, 1));
+ };
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ auto f = [](const R3& X) -> R2x2 {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return R2x2{x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1, 2 * x + 3 * y - 1, 2 * x * x + 3 * x * y - z * z,
+ 3 - 2 * x + 3 * y};
+ };
+
+ std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
+ mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
+ mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(*hybrid_mesh)));
+
+ for (auto mesh_info : mesh_list) {
+ auto mesh_name = mesh_info.first;
+ auto mesh = mesh_info.second;
+
+ Array<const R2x2> int_f_per_edge = [=] {
+ Array<R2x2> int_f(mesh->numberOfEdges());
+ auto edge_to_node_matrix = mesh->connectivity().edgeToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+ auto edge_node_list = edge_to_node_matrix[edge_id];
+ auto xr = mesh->xr();
+ R2x2 integral = zero;
+
+ LineTransformation<3> T(xr[edge_node_list[0]], xr[edge_node_list[1]]);
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.velocityNorm() * f(T(xi));
+ }
+
+ int_f[edge_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION(mesh_name)
{
- SECTION("all edges")
+ SECTION("direct formula")
{
- SECTION("EdgeValue")
+ SECTION("all edges")
{
- EdgeValue<double> values(mesh->connectivity());
- EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
- *mesh, values);
+ SECTION("EdgeValue")
+ {
+ EdgeValue<R2x2> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfEdges());
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
- EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<R2x2> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("edge list")
{
- SmallArray<double> values(mesh->numberOfEdges());
- EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ SECTION("Array")
+ {
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
- double error = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
- error += std::abs(int_f_per_edge[edge_id] - values[edge_id]);
+ REQUIRE(k == edge_list.size());
+ }
+
+ Array<R2x2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<R2x2> values = EdgeIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
- SECTION("edge list")
+ SECTION("tensorial formula")
{
- SECTION("Array")
+ SECTION("all edges")
{
- Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
-
+ SECTION("EdgeValue")
{
- size_t k = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
- edge_list[k] = edge_id;
+ EdgeValue<R2x2> values(mesh->connectivity());
+ EdgeIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
}
- REQUIRE(k == edge_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- Array<double> values =
- EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfEdges());
+
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < edge_list.size(); ++i) {
- error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<R2x2> values(mesh->numberOfEdges());
+ EdgeIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++edge_id) {
+ error += l2Norm(int_f_per_edge[edge_id] - values[edge_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("edge list")
{
- SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
-
+ SECTION("Array")
{
- size_t k = 0;
- for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
- edge_list[k] = edge_id;
+ Array<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
}
- REQUIRE(k == edge_list.size());
- }
+ Array<R2x2> values =
+ EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
- SmallArray<double> values =
- EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
- double error = 0;
- for (size_t i = 0; i < edge_list.size(); ++i) {
- error += std::abs(int_f_per_edge[edge_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<EdgeId> edge_list{mesh->numberOfEdges() / 2 + mesh->numberOfEdges() % 2};
+
+ {
+ size_t k = 0;
+ for (EdgeId edge_id = 0; edge_id < mesh->numberOfEdges(); ++(++edge_id), ++k) {
+ edge_list[k] = edge_id;
+ }
+
+ REQUIRE(k == edge_list.size());
+ }
+
+ SmallArray<R2x2> values =
+ EdgeIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, edge_list);
+
+ double error = 0;
+ for (size_t i = 0; i < edge_list.size(); ++i) {
+ error += l2Norm(int_f_per_edge[edge_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
}