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