diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
index e3a73a89ca5bcc6f5fda9039be3900805ef32843..6960901d335f894e00a86239772f33f3abc50c67 100644
--- a/src/scheme/CellIntegrator.hpp
+++ b/src/scheme/CellIntegrator.hpp
@@ -134,10 +134,12 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
@@ -161,10 +163,12 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
} else {
static_assert(Dimension == 3);
@@ -194,8 +198,10 @@ class CellIntegrator
value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * function(T(xi));
}
} else {
+ // LCOV_EXCL_START
throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
std::to_string(cell_to_node.size()) + ")");
+ // LCOV_EXCL_STOP
}
break;
}
@@ -247,8 +253,10 @@ class CellIntegrator
}
}
} else {
+ // LCOV_EXCL_START
throw NotImplementedError("integration on diamond with non-quadrangular base (" +
std::to_string(cell_to_node.size()) + ")");
+ // LCOV_EXCL_STOP
}
break;
}
@@ -275,19 +283,23 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
}
});
+ // LCOV_EXCL_START
if (invalid_cell.first) {
std::ostringstream os;
os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
throw UnexpectedError(os.str());
}
+ // LCOV_EXCL_STOP
}
template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
@@ -340,10 +352,12 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
@@ -368,10 +382,12 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
} else {
static_assert(Dimension == 3);
@@ -399,7 +415,9 @@ class CellIntegrator
value[index] += qf.weight(i_point) * T.jacobianDeterminant(xi) * function(T(xi));
}
} else {
+ // LCOV_EXCL_START
throw NotImplementedError("integration on pyramid with non-quadrangular base");
+ // LCOV_EXCL_STOP
}
break;
}
@@ -445,8 +463,10 @@ class CellIntegrator
}
}
} else {
+ // LCOV_EXCL_START
throw NotImplementedError("integration on diamond with non-quadrangular base (" +
std::to_string(cell_to_node.size()) + ")");
+ // LCOV_EXCL_STOP
}
break;
}
@@ -473,19 +493,23 @@ class CellIntegrator
}
break;
}
+ // LCOV_EXCL_START
default: {
invalid_cell = std::make_pair(true, cell_id);
break;
}
+ // LCOV_EXCL_STOP
}
}
});
+ // LCOV_EXCL_START
if (invalid_cell.first) {
std::ostringstream os;
os << "invalid cell type for integration: " << name(cell_type[invalid_cell.second]);
throw UnexpectedError(os.str());
}
+ // LCOV_EXCL_STOP
}
public:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 7f4c6a323d865a1cc9b4c2a9c1c997e5672e1a9d..8d7595e7ca2e9809829f9636e6d52de4c343910b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -53,6 +53,7 @@ add_executable (unit_tests
test_BuiltinFunctionEmbedderTable.cpp
test_BuiltinFunctionProcessor.cpp
test_CastArray.cpp
+ test_CellIntegrator.cpp
test_ConsoleManager.cpp
test_CG.cpp
test_ContinueProcessor.cpp
@@ -127,7 +128,6 @@ add_executable (unit_tests
add_executable (mpi_unit_tests
mpi_test_main.cpp
- test_CellIntegrator.cpp
test_DiscreteFunctionInterpoler.cpp
test_DiscreteFunctionP0.cpp
test_DiscreteFunctionP0Vector.cpp
diff --git a/tests/test_CellIntegrator.cpp b/tests/test_CellIntegrator.cpp
index ef0fca0b93b62710eed6bcac48a55e889824d075..97bf0b98fa9beaace5bea6bf427b69dbe75548fc 100644
--- a/tests/test_CellIntegrator.cpp
+++ b/tests/test_CellIntegrator.cpp
@@ -1,8 +1,10 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/Mesh.hpp>
#include <scheme/CellIntegrator.hpp>
@@ -480,7 +482,7 @@ TEST_CASE("CellIntegrator", "[scheme]")
{
using R3 = TinyVector<3>;
- const auto mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+ auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
auto f = [](const R3& X) -> double {
const double x = X[0];
@@ -489,261 +491,442 @@ TEST_CASE("CellIntegrator", "[scheme]")
return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
};
- Array<const double> int_f_per_cell = [=] {
- Array<double> int_f(mesh->numberOfCells());
- auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
- auto cell_type = mesh->connectivity().cellType();
-
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto cell_node_list = cell_to_node_matrix[cell_id];
- auto xr = mesh->xr();
- double integral = 0;
+ 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", DiamondDualMeshManager::instance().getDiamondDualMesh(hybrid_mesh)));
- switch (cell_type[cell_id]) {
- case CellType::Tetrahedron: {
- TetrahedronTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
- xr[cell_node_list[3]]);
- auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+ for (auto mesh_info : mesh_list) {
+ auto mesh_name = mesh_info.first;
+ auto mesh = mesh_info.second;
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi));
+ SECTION(mesh_name)
+ {
+ SECTION("direct formula")
+ {
+ Array<const double> int_f_per_cell = [=] {
+ Array<double> int_f(mesh->numberOfCells());
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ auto cell_type = mesh->connectivity().cellType();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]]);
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant() * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ PyramidTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]]);
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Prism: {
+ PrismTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]]);
+ auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
+ xr[cell_node_list[6]], xr[cell_node_list[7]]);
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4));
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_node_list.size() == 5) {
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(4));
+ { // top tetrahedron
+ TetrahedronTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+ xr[cell_node_list[4]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T0.jacobianDeterminant() * f(T0(xi));
+ }
+ }
+ { // bottom tetrahedron
+ TetrahedronTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
+ xr[cell_node_list[0]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T1.jacobianDeterminant() * f(T1(xi));
+ }
+ }
+ } else if (cell_node_list.size() == 6) {
+ auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
+ { // top pyramid
+ PyramidTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+ xr[cell_node_list[4]], xr[cell_node_list[5]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+ }
+ }
+ { // bottom pyramid
+ PyramidTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+ xr[cell_node_list[1]], xr[cell_node_list[0]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+ }
+ }
+ } else {
+ INFO("Diamond cells with more than 6 vertices are not tested");
+ REQUIRE(false);
+ }
+ break;
+ }
+ default: {
+ INFO("Diamond cells not tested yet");
+ REQUIRE(cell_type[cell_id] != CellType::Diamond);
+ }
+ }
+ int_f[cell_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("all cells")
+ {
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh,
+ values);
+
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == 0);
}
- break;
- }
- case CellType::Pyramid: {
- PyramidTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
- xr[cell_node_list[3]], xr[cell_node_list[4]]);
- auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
- }
- break;
- }
- case CellType::Prism: {
- PrismTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
- xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]]);
- auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor(4));
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
- }
- break;
- }
- case CellType::Hexahedron: {
- CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
- xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
- xr[cell_node_list[6]], xr[cell_node_list[7]]);
- auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(4));
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
- for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
- const auto& xi = qf.point(i);
- integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == 0);
}
- break;
- }
- default: {
- throw UnexpectedError("invalid cell type in 2d");
- }
- }
- int_f[cell_id] = integral;
- });
- return int_f;
- }();
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
- SECTION("direct formula")
- {
- SECTION("all cells")
- {
- SECTION("CellValue")
- {
- CellValue<double> values(mesh->connectivity());
- CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussQuadratureDescriptor(4), *mesh, values);
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ REQUIRE(error == 0);
+ }
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfCells());
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
- CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ REQUIRE(k == cell_list.size());
+ }
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
- }
+ Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
- SECTION("SmallArray")
- {
- SmallArray<double> values(mesh->numberOfCells());
- CellIntegrator::integrateTo(f, GaussQuadratureDescriptor(4), *mesh, values);
+ REQUIRE(error == 0);
+ }
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
- }
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
- }
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
- SECTION("cell list")
- {
- SECTION("Array")
- {
- Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+ REQUIRE(k == cell_list.size());
+ }
- {
- size_t k = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
- cell_list[k] = cell_id;
- }
+ SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
- REQUIRE(k == cell_list.size());
- }
-
- Array<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
- double error = 0;
- for (size_t i = 0; i < cell_list.size(); ++i) {
- error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ REQUIRE(error == 0);
+ }
}
-
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- SECTION("SmallArray")
+ SECTION("tensorial formula")
{
- SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
-
+ Array<const double> int_f_per_cell = [=] {
+ Array<double> int_f(mesh->numberOfCells());
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ auto cell_type = mesh->connectivity().cellType();
+
+ auto qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ auto xr = mesh->xr();
+ double integral = 0;
+
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[3]],
+ xr[cell_node_list[3]], xr[cell_node_list[3]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Pyramid: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
+ xr[cell_node_list[4]], xr[cell_node_list[4]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Prism: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[2]], xr[cell_node_list[3]], xr[cell_node_list[4]],
+ xr[cell_node_list[5]], xr[cell_node_list[5]]);
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Hexahedron: {
+ CubeTransformation T(xr[cell_node_list[0]], xr[cell_node_list[1]], xr[cell_node_list[2]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[5]],
+ xr[cell_node_list[6]], xr[cell_node_list[7]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T.jacobianDeterminant(xi) * f(T(xi));
+ }
+ break;
+ }
+ case CellType::Diamond: {
+ if (cell_node_list.size() == 5) {
+ { // top tetrahedron
+ CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+ xr[cell_node_list[3]], xr[cell_node_list[4]], xr[cell_node_list[4]],
+ xr[cell_node_list[4]], xr[cell_node_list[4]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+ }
+ }
+ { // bottom tetrahedron
+ CubeTransformation T1(xr[cell_node_list[3]], xr[cell_node_list[2]], xr[cell_node_list[1]],
+ xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+ xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+ }
+ }
+ } else if (cell_node_list.size() == 6) {
+ { // top pyramid
+ CubeTransformation T0(xr[cell_node_list[1]], xr[cell_node_list[2]], xr[cell_node_list[3]],
+ xr[cell_node_list[4]], xr[cell_node_list[5]], xr[cell_node_list[5]],
+ xr[cell_node_list[5]], xr[cell_node_list[5]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T0.jacobianDeterminant(xi) * f(T0(xi));
+ }
+ }
+ { // bottom pyramid
+ CubeTransformation T1(xr[cell_node_list[4]], xr[cell_node_list[3]], xr[cell_node_list[2]],
+ xr[cell_node_list[1]], xr[cell_node_list[0]], xr[cell_node_list[0]],
+ xr[cell_node_list[0]], xr[cell_node_list[0]]);
+
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const auto& xi = qf.point(i);
+ integral += qf.weight(i) * T1.jacobianDeterminant(xi) * f(T1(xi));
+ }
+ }
+ } else {
+ INFO("Diamond cells with more than 6 vertices are not tested");
+ REQUIRE(false);
+ }
+ break;
+ }
+ default: {
+ INFO("Diamond cells not tested yet");
+ REQUIRE(cell_type[cell_id] != CellType::Diamond);
+ }
+ }
+ int_f[cell_id] = integral;
+ });
+
+ return int_f;
+ }();
+
+ SECTION("all cells")
{
- size_t k = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
- cell_list[k] = cell_id;
+ SECTION("CellValue")
+ {
+ CellValue<double> values(mesh->connectivity());
+ CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLegendreQuadratureDescriptor(10),
+ *mesh, values);
+
+ auto cell_type = mesh->connectivity().cellType();
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
+
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(k == cell_list.size());
- }
-
- SmallArray<double> values = CellIntegrator::integrate(f, GaussQuadratureDescriptor(4), *mesh, cell_list);
+ SECTION("Array")
+ {
+ Array<double> values(mesh->numberOfCells());
- double error = 0;
- for (size_t i = 0; i < cell_list.size(); ++i) {
- error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
- }
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
- }
- }
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
- SECTION("tensorial formula")
- {
- SECTION("all cells")
- {
- SECTION("CellValue")
- {
- CellValue<double> values(mesh->connectivity());
- CellIntegrator::integrateTo([=](const R3 x) { return f(x); }, GaussLobattoQuadratureDescriptor(4), *mesh,
- values);
-
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
- }
-
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
- SECTION("Array")
- {
- Array<double> values(mesh->numberOfCells());
+ SECTION("SmallArray")
+ {
+ SmallArray<double> values(mesh->numberOfCells());
+ CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
- CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ double error = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ }
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
+ }
}
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SECTION("cell list")
+ {
+ SECTION("Array")
+ {
+ Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
- SECTION("SmallArray")
- {
- SmallArray<double> values(mesh->numberOfCells());
- CellIntegrator::integrateTo(f, GaussLobattoQuadratureDescriptor(4), *mesh, values);
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
- double error = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- error += std::abs(int_f_per_cell[cell_id] - values[cell_id]);
- }
+ REQUIRE(k == cell_list.size());
+ }
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
- }
+ Array<double> values =
+ CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
- SECTION("cell list")
- {
- SECTION("Array")
- {
- Array<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
- {
- size_t k = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
- cell_list[k] = cell_id;
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
- REQUIRE(k == cell_list.size());
- }
+ SECTION("SmallArray")
+ {
+ SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
- Array<double> values = CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
+ {
+ size_t k = 0;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
+ cell_list[k] = cell_id;
+ }
- double error = 0;
- for (size_t i = 0; i < cell_list.size(); ++i) {
- error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
- }
+ REQUIRE(k == cell_list.size());
+ }
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
- }
+ SmallArray<double> values =
+ CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
- SECTION("SmallArray")
- {
- SmallArray<CellId> cell_list{mesh->numberOfCells() / 2 + mesh->numberOfCells() % 2};
+ double error = 0;
+ for (size_t i = 0; i < cell_list.size(); ++i) {
+ error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
+ }
- {
- size_t k = 0;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++(++cell_id), ++k) {
- cell_list[k] = cell_id;
+ REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
-
- REQUIRE(k == cell_list.size());
}
-
- SmallArray<double> values =
- CellIntegrator::integrate(f, GaussLobattoQuadratureDescriptor(4), *mesh, cell_list);
-
- double error = 0;
- for (size_t i = 0; i < cell_list.size(); ++i) {
- error += std::abs(int_f_per_cell[cell_list[i]] - values[i]);
- }
-
- REQUIRE(error == Catch::Approx(0).margin(1E-10));
}
}
}
diff --git a/tests/test_CubeTransformation.cpp b/tests/test_CubeTransformation.cpp
index 2de108f9c96c5c3659c20b4c3e7d314d5342280c..7d060ccd6faf409fb08813607209da1d9da616b1 100644
--- a/tests/test_CubeTransformation.cpp
+++ b/tests/test_CubeTransformation.cpp
@@ -2,6 +2,8 @@
#include <catch2/catch_test_macros.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
@@ -11,88 +13,88 @@ TEST_CASE("CubeTransformation", "[geometry]")
{
using R3 = TinyVector<3>;
- const R3 a_hat = {-1, -1, -1};
- const R3 b_hat = {+1, -1, -1};
- const R3 c_hat = {+1, +1, -1};
- const R3 d_hat = {-1, +1, -1};
- const R3 e_hat = {-1, -1, +1};
- const R3 f_hat = {+1, -1, +1};
- const R3 g_hat = {+1, +1, +1};
- const R3 h_hat = {-1, +1, +1};
-
- const R3 m_hat = zero;
+ SECTION("invertible transformation")
+ {
+ const R3 a_hat = {-1, -1, -1};
+ const R3 b_hat = {+1, -1, -1};
+ const R3 c_hat = {+1, +1, -1};
+ const R3 d_hat = {-1, +1, -1};
+ const R3 e_hat = {-1, -1, +1};
+ const R3 f_hat = {+1, -1, +1};
+ const R3 g_hat = {+1, +1, +1};
+ const R3 h_hat = {-1, +1, +1};
- const R3 a = {0, 0, 0};
- const R3 b = {3, 1, 3};
- const R3 c = {2, 5, 2};
- const R3 d = {0, 3, 1};
- const R3 e = {1, 2, 5};
- const R3 f = {3, 1, 7};
- const R3 g = {2, 5, 5};
- const R3 h = {0, 3, 6};
+ const R3 m_hat = zero;
- const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
+ const R3 a = {0, 0, 0};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {0, 3, 1};
+ const R3 e = {1, 2, 5};
+ const R3 f = {3, 1, 7};
+ const R3 g = {2, 5, 5};
+ const R3 h = {0, 3, 6};
- const CubeTransformation t(a, b, c, d, e, f, g, h);
+ const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
- SECTION("values")
- {
- REQUIRE(l2Norm(t(a_hat) - a) == Catch::Approx(0));
- REQUIRE(l2Norm(t(b_hat) - b) == Catch::Approx(0));
- REQUIRE(l2Norm(t(c_hat) - c) == Catch::Approx(0));
- REQUIRE(l2Norm(t(d_hat) - d) == Catch::Approx(0));
- REQUIRE(l2Norm(t(e_hat) - e) == Catch::Approx(0));
- REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
- REQUIRE(l2Norm(t(g_hat) - g) == Catch::Approx(0));
- REQUIRE(l2Norm(t(h_hat) - h) == Catch::Approx(0));
-
- REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0));
- }
+ const CubeTransformation t(a, b, c, d, e, f, g, h);
- SECTION("Jacobian determinant")
- {
- SECTION("at points")
+ SECTION("values")
{
- auto detJ = [](const R3& X) {
- const double x = X[0];
- const double y = X[1];
- const double z = X[2];
-
- return -(3 * x * y * y * z + 9 * y * y * z - 6 * x * x * y * z + 6 * x * y * z - 96 * y * z + 14 * x * x * z -
- 49 * x * z + 119 * z - 21 * x * y * y + 21 * y * y + 90 * x * y - 10 * y + 40 * x * x - 109 * x -
- 491) /
- 128;
- };
+ REQUIRE(l2Norm(t(a_hat) - a) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(b_hat) - b) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(c_hat) - c) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(d_hat) - d) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(e_hat) - e) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(g_hat) - g) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(h_hat) - h) == Catch::Approx(0));
- REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
- REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
- REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
- REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
- REQUIRE(t.jacobianDeterminant(e_hat) == Catch::Approx(detJ(e_hat)));
- REQUIRE(t.jacobianDeterminant(f_hat) == Catch::Approx(detJ(f_hat)));
- REQUIRE(t.jacobianDeterminant(g_hat) == Catch::Approx(detJ(g_hat)));
- REQUIRE(t.jacobianDeterminant(h_hat) == Catch::Approx(detJ(h_hat)));
-
- REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+ REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0));
}
- SECTION("Gauss order 3")
+ SECTION("Jacobian determinant")
{
- // The jacobian determinant is of maximal degree 2 in each variable
- const QuadratureFormula<3>& gauss =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(2));
+ SECTION("at points")
+ {
+ auto detJ = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
- double volume = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ return -(3 * x * y * y * z + 9 * y * y * z - 6 * x * x * y * z + 6 * x * y * z - 96 * y * z + 14 * x * x * z -
+ 49 * x * z + 119 * z - 21 * x * y * y + 21 * y * y + 90 * x * y - 10 * y + 40 * x * x - 109 * x -
+ 491) /
+ 128;
+ };
+
+ REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
+ REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
+ REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
+ REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
+ REQUIRE(t.jacobianDeterminant(e_hat) == Catch::Approx(detJ(e_hat)));
+ REQUIRE(t.jacobianDeterminant(f_hat) == Catch::Approx(detJ(f_hat)));
+ REQUIRE(t.jacobianDeterminant(g_hat) == Catch::Approx(detJ(g_hat)));
+ REQUIRE(t.jacobianDeterminant(h_hat) == Catch::Approx(detJ(h_hat)));
+
+ REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
}
- // 353. / 12 is actually the volume of the hexahedron
- REQUIRE(volume == Catch::Approx(353. / 12));
- }
+ SECTION("Gauss order 3")
+ {
+ // The jacobian determinant is of maximal degree 2 in each variable
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(2));
+
+ double volume = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ }
+
+ // 353. / 12 is actually the volume of the hexahedron
+ REQUIRE(volume == Catch::Approx(353. / 12));
+ }
- SECTION("Gauss order 4")
- {
auto p = [](const R3& X) {
const double& x = X[0];
const double& y = X[1];
@@ -101,17 +103,196 @@ TEST_CASE("CubeTransformation", "[geometry]")
return 2 * x * x + 3 * y * z + z * z + 2 * x - 3 * y + 0.5 * z + 3;
};
- // Jacbian determinant is a degree 2 polynomial, so the
- // following formula is required to reach exactness
+ SECTION("Gauss order 6")
+ {
+ // Jacbian determinant is a degree 2 polynomial, so the
+ // following formula is required to reach exactness
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(6));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(488879. / 360));
+ }
+
+ SECTION("Gauss-Legendre order 4")
+ {
+ // Jacbian determinant is a degree 2 polynomial, so the
+ // following formula is required to reach exactness
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(488879. / 360));
+ }
+
+ SECTION("Gauss-Lobatto order 4")
+ {
+ // Jacbian determinant is a degree 2 polynomial, so the
+ // following formula is required to reach exactness
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(488879. / 360));
+ }
+ }
+ }
+
+ SECTION("degenerate to tetrahedron")
+ {
+ using R3 = TinyVector<3>;
+
+ const R3 a = {1, 2, 1};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {2, 3, 4};
+
+ const CubeTransformation t(a, b, c, c, d, d, d, d);
+
+ auto p = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return 2 * x * x + 3 * x * y + y * y + 3 * y - z * z + 2 * x * z + 2 * z;
+ };
+
+ SECTION("Polynomial Gauss integral")
+ {
+ QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(6));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const R3 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(231. / 2));
+ }
+
+ SECTION("Polynomial Gauss-Lobatto integral")
+ {
+ QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const R3 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(231. / 2));
+ }
+
+ SECTION("Polynomial Gauss-Legendre integral")
+ {
+ QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ const R3 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(231. / 2));
+ }
+ }
+
+ SECTION("degenerate to prism")
+ {
+ const R3 a = {1, 2, 0};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {0, 3, 1};
+ const R3 e = {1, 2, 5};
+ const R3 f = {3, 1, 7};
+
+ const CubeTransformation t(a, b, c, c, d, e, f, f);
+
+ auto p = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+
+ return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
+ };
+
+ SECTION("Polynomial Gauss integral")
+ {
+ // 8 is the minimum quadrature rule to integrate the polynomial
+ // on the prism due to the jacobian determinant expression
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(8));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(30377. / 90));
+ }
+
+ SECTION("Polynomial Gauss-Legendre integral")
+ {
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
- REQUIRE(integral == Catch::Approx(8795. / 72));
+ REQUIRE(integral == Catch::Approx(30377. / 90));
}
+
+ SECTION("Polynomial Gauss-Lobatto integral")
+ {
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(30377. / 90));
+ }
+ }
+
+ SECTION("degenerate to pyramid")
+ {
+ const R3 a = {1, 2, 0};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {0, 3, 1};
+ const R3 e = {1, 2, 5};
+
+ const CubeTransformation t(a, b, c, d, e, e, e, e);
+
+ auto p = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+
+ return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
+ };
+
+ // 4 is the minimum quadrature rule to integrate the polynomial on the pyramid
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(20));
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
+ }
+
+ REQUIRE(integral == Catch::Approx(7238. / 15));
}
}
diff --git a/tests/test_LineTransformation.cpp b/tests/test_LineTransformation.cpp
index 764a842b8fdc0d817955d4d0e4e1c7ef4b5e2e48..908d2d693973d13d48288cda7f7d021ee187c7a7 100644
--- a/tests/test_LineTransformation.cpp
+++ b/tests/test_LineTransformation.cpp
@@ -23,6 +23,25 @@ TEST_CASE("LineTransformation", "[geometry]")
REQUIRE(t(1)[0] == Catch::Approx(2.7));
REQUIRE(t.jacobianDeterminant() == Catch::Approx(1.2));
+
+ auto p = [](const R1& X) {
+ const double x = X[0];
+ return 2 * x * x + 3 * x + 2;
+ };
+
+ QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ sum += qf.weight(i) * t.jacobianDeterminant() * p(t(qf.point(i)));
+ }
+
+ auto P = [](const R1& X) {
+ const double x = X[0];
+ return 2. / 3 * x * x * x + 3. / 2 * x * x + 2 * x;
+ };
+
+ REQUIRE(sum == Catch::Approx(P(b) - P(a)));
}
SECTION("2D")
diff --git a/tests/test_PrismTransformation.cpp b/tests/test_PrismTransformation.cpp
index e170da3c05cf2643ab22958509b2cf62039f3fa1..f933fe308fbc64915faf04449aeca364f0904053 100644
--- a/tests/test_PrismTransformation.cpp
+++ b/tests/test_PrismTransformation.cpp
@@ -97,10 +97,10 @@ TEST_CASE("PrismTransformation", "[geometry]")
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
- REQUIRE(integral == Catch::Approx(17627. / 720));
+ REQUIRE(integral == Catch::Approx(30377. / 90));
}
}
}
diff --git a/tests/test_PyramidTransformation.cpp b/tests/test_PyramidTransformation.cpp
index ca51bf080045e0df06d706b60bc5fd1dca497d44..23b7369e625b20307ddab31d222b1b2ae3d607b4 100644
--- a/tests/test_PyramidTransformation.cpp
+++ b/tests/test_PyramidTransformation.cpp
@@ -1,8 +1,10 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
#include <geometry/PyramidTransformation.hpp>
// clazy:excludeall=non-pod-global-static
@@ -83,14 +85,14 @@ TEST_CASE("PyramidTransformation", "[geometry]")
return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
};
- // 3 is the minimum quadrature rule to integrate the polynomial on the pyramid
- const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(3));
+ // 4 is the minimum quadrature rule to integrate the polynomial on the pyramid
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
- REQUIRE(integral == Catch::Approx(1513. / 60));
+ REQUIRE(integral == Catch::Approx(213095. / 448));
}
}
}
diff --git a/tests/test_SquareTransformation.cpp b/tests/test_SquareTransformation.cpp
index fa12d756f27945d0189f19785c7c5bc1caccb0d5..a1ead02fc75d25f5ff0194b4157e22497d6005dc 100644
--- a/tests/test_SquareTransformation.cpp
+++ b/tests/test_SquareTransformation.cpp
@@ -2,6 +2,8 @@
#include <catch2/catch_test_macros.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/SquareTransformation.hpp>
@@ -97,10 +99,69 @@ TEST_CASE("SquareTransformation", "[geometry]")
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
- REQUIRE(integral == Catch::Approx(-479. / 2));
+ REQUIRE(integral == Catch::Approx(-76277. / 24));
+ }
+ }
+ }
+
+ SECTION("degenerate 2D ")
+ {
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ const R2 a = {1, 2};
+ const R2 b = {3, 1};
+ const R2 c = {2, 5};
+
+ const SquareTransformation<2> t(a, b, c, c);
+
+ auto p = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return 2 * x * x + 3 * x * y + y * y + 3 * y + 1;
+ };
+
+ SECTION("Gauss")
+ {
+ QuadratureFormula<2> qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ R2 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(3437. / 24));
+ }
+
+ SECTION("Gauss Legendre")
+ {
+ QuadratureFormula<2> qf = QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ R2 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(3437. / 24));
+ }
+
+ SECTION("Gauss Lobatto")
+ {
+ QuadratureFormula<2> qf = QuadratureManager::instance().getSquareFormula(GaussLobattoQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ R2 xi = qf.point(i);
+ sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
+ }
+
+ REQUIRE(sum == Catch::Approx(3437. / 24));
}
}
}
diff --git a/tests/test_TetrahedronTransformation.cpp b/tests/test_TetrahedronTransformation.cpp
index 1a9d20275fd1df42db4e5ea723b585c2872e69c4..91739b99e7030e1b8490869fb770f1978cce9334 100644
--- a/tests/test_TetrahedronTransformation.cpp
+++ b/tests/test_TetrahedronTransformation.cpp
@@ -1,6 +1,8 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/TetrahedronTransformation.hpp>
// clazy:excludeall=non-pod-global-static
@@ -37,4 +39,23 @@ TEST_CASE("TetrahedronTransformation", "[geometry]")
REQUIRE(t({0.25, 0.25, 0.25})[2] == Catch::Approx(2.5));
REQUIRE(t.jacobianDeterminant() == Catch::Approx(14));
+
+ SECTION("Polynomial integral")
+ {
+ auto p = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return 2 * x * x + 3 * x * y + y * y + 3 * y - z * z + 2 * x * z + 2 * z;
+ };
+
+ QuadratureFormula<3> qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ sum += qf.weight(i) * t.jacobianDeterminant() * p(t(qf.point(i)));
+ }
+
+ REQUIRE(sum == Catch::Approx(231. / 2));
+ }
}
diff --git a/tests/test_TriangleTransformation.cpp b/tests/test_TriangleTransformation.cpp
index 99330e41d7376d7a7265ba19798d6a0af82c3faa..e9561e859d727ea6c4019eaa78de125fbab8ba00 100644
--- a/tests/test_TriangleTransformation.cpp
+++ b/tests/test_TriangleTransformation.cpp
@@ -32,6 +32,24 @@ TEST_CASE("TriangleTransformation", "[geometry]")
REQUIRE(t({1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
REQUIRE(t.jacobianDeterminant() == Catch::Approx(7));
+
+ SECTION("Polynomial integral")
+ {
+ auto p = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return 2 * x * x + 3 * x * y + y * y + 3 * y + 1;
+ };
+
+ QuadratureFormula<2> qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(2));
+
+ double sum = 0;
+ for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+ sum += qf.weight(i) * t.jacobianDeterminant() * p(t(qf.point(i)));
+ }
+
+ REQUIRE(sum == Catch::Approx(3437. / 24));
+ }
}
SECTION("3D")