diff --git a/src/scheme/FaceIntegrator.hpp b/src/scheme/FaceIntegrator.hpp
index 30ee39077ea12d898b4bdee00f12970be79e8753..657ca7a88a3a959f2856bc34ae8938ab91561785 100644
--- a/src/scheme/FaceIntegrator.hpp
+++ b/src/scheme/FaceIntegrator.hpp
@@ -144,9 +144,11 @@ class FaceIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_face = std::make_pair(true, face_id);
}
+ // LCOV_EXCL_STOP
}
});
@@ -159,7 +161,9 @@ class FaceIntegrator
// LCOV_EXCL_STOP
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("invalid dimension for face integration");
+ // LCOV_EXCL_STOP
}
}
diff --git a/tests/test_FaceIntegrator.cpp b/tests/test_FaceIntegrator.cpp
index 184f1a8f97db137f99f638f5feef048f4d78ec04..945d765a6d0a530c7c7babb5d8b5800d5af033d3 100644
--- a/tests/test_FaceIntegrator.cpp
+++ b/tests/test_FaceIntegrator.cpp
@@ -15,487 +15,1477 @@
TEST_CASE("FaceIntegrator", "[scheme]")
{
- SECTION("2D")
+ SECTION("scalar")
{
- using R2 = TinyVector<2>;
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
- const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+ const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
- auto f = [](const R2& X) -> double {
- const double x = X[0];
- const double y = X[1];
- return x * x + 2 * x * y + 3 * y * y + 2;
- };
+ auto f = [](const R2& X) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ return x * x + 2 * x * y + 3 * y * y + 2;
+ };
- Array<const double> int_f_per_face = [=] {
- Array<double> int_f(mesh->numberOfFaces());
- auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ Array<const double> int_f_per_face = [=] {
+ Array<double> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
- parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
- auto face_node_list = face_to_node_matrix[face_id];
- auto xr = mesh->xr();
- double integral = 0;
- LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]);
- auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+ LineTransformation<2> T(xr[face_node_list[0]], xr[face_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));
- }
+ 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[face_id] = integral;
- });
+ int_f[face_id] = integral;
+ });
- return int_f;
- }();
+ return int_f;
+ }();
- SECTION("direct formula")
- {
- SECTION("all faces")
+ SECTION("direct formula")
{
- SECTION("FaceValue")
+ SECTION("all faces")
{
- FaceValue<double> values(mesh->connectivity());
- FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
- FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<double> values(mesh->numberOfFaces());
- FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
- SECTION("face list")
+ SECTION("tensorial formula")
{
- SECTION("Array")
+ SECTION("all faces")
{
- Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
-
+ SECTION("FaceValue")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
}
- REQUIRE(k == face_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
}
- REQUIRE(k == face_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
+ }
+ }
+ }
+
+ 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_face = [=] {
+ Array<double> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (face_node_list.size()) {
+ case 3: {
+ TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]);
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]],
+ xr[face_node_list[3]]);
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid face (node number must be 3 or 4)");
+ }
+ }
+ int_f[face_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
- SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values =
+ FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("tensorial formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<double> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<double> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
}
}
}
+ }
- SECTION("tensorial formula")
+ SECTION("R^d")
+ {
+ SECTION("2D")
{
- SECTION("all faces")
+ using R2 = TinyVector<2>;
+
+ const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ auto f = [](const R2& X) -> R2 {
+ const double x = X[0];
+ const double y = X[1];
+ return R2{x * x + 2 * x * y + 3 * y * y + 2, 3 * x - 2 * y};
+ };
+
+ Array<const R2> int_f_per_face = [=] {
+ Array<R2> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ R2 integral = zero;
+ LineTransformation<2> T(xr[face_node_list[0]], xr[face_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[face_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("direct formula")
{
- SECTION("FaceValue")
+ SECTION("all faces")
{
- FaceValue<double> values(mesh->connectivity());
- FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
- values);
+ SECTION("FaceValue")
+ {
+ FaceValue<R2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2& x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("Array")
+ {
+ Array<R2> values(mesh->numberOfFaces());
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
- FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<R2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<double> values(mesh->numberOfFaces());
- FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
- SECTION("face list")
+ SECTION("tensorial formula")
{
- SECTION("Array")
+ SECTION("all faces")
{
- Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
-
+ SECTION("FaceValue")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ FaceValue<R2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
}
- REQUIRE(k == face_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- Array<double> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+ SECTION("Array")
+ {
+ Array<R2> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<R2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+ Array<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
}
- REQUIRE(k == face_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ using R2 = TinyVector<2>;
+ 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, 3 * x - 2 * y + z};
+ };
+
+ 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_face = [=] {
+ Array<R2> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ R2 integral = zero;
+
+ switch (face_node_list.size()) {
+ case 3: {
+ TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]);
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]],
+ xr[face_node_list[3]]);
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid face (node number must be 3 or 4)");
+ }
+ }
+ int_f[face_id] = integral;
+ });
- SmallArray<double> values =
- FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+ return int_f;
+ }();
+
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<R2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<R2> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ SECTION("SmallArray")
+ {
+ SmallArray<R2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("tensorial formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
+ {
+ FaceValue<R2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("Array")
+ {
+ Array<R2> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<R2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<R2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+ }
}
}
}
}
- SECTION("3D")
+ SECTION("R^dxd")
{
- 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;
+ 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));
};
- 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)));
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ const auto mesh = MeshDataBaseForTests::get().hybrid2DMesh();
- for (auto mesh_info : mesh_list) {
- auto mesh_name = mesh_info.first;
- auto mesh = mesh_info.second;
+ auto f = [](const R2& X) -> R2x2 {
+ const double x = X[0];
+ const double y = X[1];
+ return R2x2{x * x + 2 * x * y + 3 * y * y + 2, 3 * x - 2 * y, 2 * x * x - 2 * y, 2 + x * y};
+ };
- Array<const double> int_f_per_face = [=] {
- Array<double> int_f(mesh->numberOfFaces());
+ Array<const R2x2> int_f_per_face = [=] {
+ Array<R2x2> int_f(mesh->numberOfFaces());
auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
parallel_for(
mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
auto face_node_list = face_to_node_matrix[face_id];
auto xr = mesh->xr();
- double integral = 0;
+ R2x2 integral = zero;
+ LineTransformation<2> T(xr[face_node_list[0]], xr[face_node_list[1]]);
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
- switch (face_node_list.size()) {
- case 3: {
- TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]);
- auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi));
- }
- break;
- }
- case 4: {
- SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]],
- xr[face_node_list[3]]);
- auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi));
- }
- break;
- }
- default: {
- throw UnexpectedError("invalid face (node number must be 3 or 4)");
- }
+ 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[face_id] = integral;
});
return int_f;
}();
- SECTION(mesh_name)
+ SECTION("direct formula")
{
- SECTION("direct formula")
+ SECTION("all faces")
{
- SECTION("all faces")
+ SECTION("FaceValue")
{
- SECTION("FaceValue")
- {
- FaceValue<double> values(mesh->connectivity());
- FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
- values);
+ FaceValue<R2x2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2& x) { return f(x); }, GaussQuadratureDescriptor(2), *mesh, values);
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
- }
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
}
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfFaces());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+
+ SECTION("SmallArray")
+ {
+ SmallArray<R2x2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
- FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ REQUIRE(k == face_list.size());
}
- SECTION("SmallArray")
- {
- SmallArray<double> values(mesh->numberOfFaces());
- FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ Array<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
}
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
+ }
+ }
- SECTION("face list")
+ SECTION("tensorial formula")
+ {
+ SECTION("all faces")
+ {
+ SECTION("FaceValue")
{
- SECTION("Array")
- {
- Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+ FaceValue<R2x2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R2 x) { return f(x); }, GaussLobattoQuadratureDescriptor(2), *mesh,
+ values);
- {
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
- }
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- REQUIRE(k == face_list.size());
- }
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- Array<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfFaces());
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
- }
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
}
- SECTION("SmallArray")
- {
- SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- {
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
- }
+ SECTION("SmallArray")
+ {
+ SmallArray<R2x2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(2), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
+ }
+
+ SECTION("face list")
+ {
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
- REQUIRE(k == face_list.size());
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
}
- SmallArray<double> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<R2x2> values = FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ REQUIRE(k == face_list.size());
}
+
+ SmallArray<R2x2> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(2), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
}
+ }
+ }
+
+ 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, 3 * x - 2 * y + z, 2 * x - y * z, 3 * z - x * x};
+ };
+
+ 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_face = [=] {
+ Array<R2x2> int_f(mesh->numberOfFaces());
+ auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ auto xr = mesh->xr();
+ R2x2 integral = zero;
+
+ switch (face_node_list.size()) {
+ case 3: {
+ TriangleTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]]);
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm() * f(T(xi));
+ }
+ break;
+ }
+ case 4: {
+ SquareTransformation<3> T(xr[face_node_list[0]], xr[face_node_list[1]], xr[face_node_list[2]],
+ xr[face_node_list[3]]);
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.areaVariationNorm(xi) * f(T(xi));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid face (node number must be 3 or 4)");
+ }
+ }
+ int_f[face_id] = integral;
+ });
+
+ return int_f;
+ }();
- SECTION("tensorial formula")
+ SECTION(mesh_name)
{
- SECTION("all faces")
+ SECTION("direct formula")
{
- SECTION("FaceValue")
+ SECTION("all faces")
{
- FaceValue<double> values(mesh->connectivity());
- FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
- *mesh, values);
+ SECTION("FaceValue")
+ {
+ FaceValue<R2x2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfFaces());
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
- FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<R2x2> values(mesh->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<double> values(mesh->numberOfFaces());
- FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ SECTION("Array")
+ {
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
- double error = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- error += std::abs(int_f_per_face[face_id] - values[face_id]);
+ REQUIRE(k == face_list.size());
+ }
+
+ Array<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2x2> values = FaceIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
- SECTION("face list")
+ SECTION("tensorial formula")
{
- SECTION("Array")
+ SECTION("all faces")
{
- Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
-
+ SECTION("FaceValue")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ FaceValue<R2x2> values(mesh->connectivity());
+ FaceIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
}
- REQUIRE(k == face_list.size());
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- Array<double> values =
- FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+ SECTION("Array")
+ {
+ Array<R2x2> values(mesh->numberOfFaces());
+
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_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->numberOfFaces());
+ FaceIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+
+ double error = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ error += l2Norm(int_f_per_face[face_id] - values[face_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- SECTION("SmallArray")
+ SECTION("face list")
{
- SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
-
+ SECTION("Array")
{
- size_t k = 0;
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
- face_list[k] = face_id;
+ Array<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
}
- REQUIRE(k == face_list.size());
- }
+ Array<R2x2> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
- SmallArray<double> values =
- FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
- double error = 0;
- for (size_t i = 0; i < face_list.size(); ++i) {
- error += std::abs(int_f_per_face[face_list[i]] - values[i]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ SECTION("SmallArray")
+ {
+ SmallArray<FaceId> face_list{mesh->numberOfFaces() / 2 + mesh->numberOfFaces() % 2};
+
+ {
+ size_t k = 0;
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++(++face_id), ++k) {
+ face_list[k] = face_id;
+ }
+
+ REQUIRE(k == face_list.size());
+ }
+
+ SmallArray<R2x2> values =
+ FaceIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, face_list);
+
+ double error = 0;
+ for (size_t i = 0; i < face_list.size(); ++i) {
+ error += l2Norm(int_f_per_face[face_list[i]] - values[i]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
}
}