From 60ecad1d7c49befa8b9aa654b0787c86b2ba0d08 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 27 Oct 2021 19:04:00 +0200
Subject: [PATCH] Change reference triangle to the classical finite element P1
one
---
src/analysis/SimplicialGaussQuadrature.cpp | 10 +++++-----
src/geometry/AffineTransformation.hpp | 8 ++++----
tests/test_AffineTransformation.cpp | 18 +++++++++---------
tests/test_EconomicalGaussQuadrature.cpp | 20 +++++++++++++++-----
tests/test_SimplicialGaussQuadrature.cpp | 14 +++++++++++---
5 files changed, 44 insertions(+), 26 deletions(-)
diff --git a/src/analysis/SimplicialGaussQuadrature.cpp b/src/analysis/SimplicialGaussQuadrature.cpp
index 5aed8d461..ab0ba3549 100644
--- a/src/analysis/SimplicialGaussQuadrature.cpp
+++ b/src/analysis/SimplicialGaussQuadrature.cpp
@@ -22,15 +22,15 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
auto fill_quadrature_points = [](auto quadrature_group_list, auto& point_list, auto& weight_list) {
Assert(point_list.size() == weight_list.size());
- const R2 A = {-1, -1};
- const R2 B = {+1, -1};
- const R2 C = {-1, +1};
+ const R2 A = {0, 0};
+ const R2 B = {1, 0};
+ const R2 C = {0, 1};
size_t k = 0;
for (size_t i = 0; i < quadrature_group_list.size(); ++i) {
- const auto [half_weight, lambda_1, lambda_2, lambda_3] = quadrature_group_list[i];
+ const auto [unit_weight, lambda_1, lambda_2, lambda_3] = quadrature_group_list[i];
- const double w = 2 * half_weight;
+ const double w = 0.5 * unit_weight;
if (lambda_2 == lambda_3) {
if (lambda_1 == lambda_2) {
diff --git a/src/geometry/AffineTransformation.hpp b/src/geometry/AffineTransformation.hpp
index e7179ef5d..2abaff4fa 100644
--- a/src/geometry/AffineTransformation.hpp
+++ b/src/geometry/AffineTransformation.hpp
@@ -47,11 +47,11 @@ class AffineTransformation
const TinyVector<Dimension>& c = image_nodes[2];
for (size_t i = 0; i < Dimension; ++i) {
- m_jacobian(i, 0) = 0.5 * (b[i] - a[i]);
- m_jacobian(i, 1) = 0.5 * (c[i] - a[i]);
-
- m_shift[i] = 0.5 * (b[i] + c[i]);
+ m_jacobian(i, 0) = b[i] - a[i];
+ m_jacobian(i, 1) = c[i] - a[i];
}
+
+ m_shift = a;
} else {
static_assert(Dimension == 3);
diff --git a/tests/test_AffineTransformation.cpp b/tests/test_AffineTransformation.cpp
index cce135741..52d962806 100644
--- a/tests/test_AffineTransformation.cpp
+++ b/tests/test_AffineTransformation.cpp
@@ -36,19 +36,19 @@ TEST_CASE("AffineTransformation", "[geometry]")
const std::array points = {a, b, c};
const AffineTransformation<2> t(points);
- REQUIRE(t({-1, -1})[0] == Catch::Approx(1));
- REQUIRE(t({-1, -1})[1] == Catch::Approx(2));
+ REQUIRE(t({0, 0})[0] == Catch::Approx(1));
+ REQUIRE(t({0, 0})[1] == Catch::Approx(2));
- REQUIRE(t({+1, -1})[0] == Catch::Approx(3));
- REQUIRE(t({+1, -1})[1] == Catch::Approx(1));
+ REQUIRE(t({1, 0})[0] == Catch::Approx(3));
+ REQUIRE(t({1, 0})[1] == Catch::Approx(1));
- REQUIRE(t({-1, +1})[0] == Catch::Approx(2));
- REQUIRE(t({-1, +1})[1] == Catch::Approx(5));
+ REQUIRE(t({0, 1})[0] == Catch::Approx(2));
+ REQUIRE(t({0, 1})[1] == Catch::Approx(5));
- REQUIRE(t({-1. / 3, -1. / 3})[0] == Catch::Approx(2));
- REQUIRE(t({-1. / 3, -1. / 3})[1] == Catch::Approx(8. / 3));
+ REQUIRE(t({1. / 3, 1. / 3})[0] == Catch::Approx(2));
+ REQUIRE(t({1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
- REQUIRE(t.jacobianDeterminant() == Catch::Approx(7. / 4));
+ REQUIRE(t.jacobianDeterminant() == Catch::Approx(7));
}
SECTION("3D")
diff --git a/tests/test_EconomicalGaussQuadrature.cpp b/tests/test_EconomicalGaussQuadrature.cpp
index 1bbe8ee7a..af86b37ff 100644
--- a/tests/test_EconomicalGaussQuadrature.cpp
+++ b/tests/test_EconomicalGaussQuadrature.cpp
@@ -2,8 +2,9 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
-#include <analysis/EconomicalGaussQuadrature.hpp>
#include <geometry/AffineTransformation.hpp>
+
+#include <analysis/EconomicalGaussQuadrature.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
@@ -538,17 +539,26 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
};
auto integrate_on_rectangle = [](auto f, auto quadrature_formula, const std::array<TinyVector<2>, 3>& triangle) {
- AffineTransformation<2> t{triangle};
+ const auto& A = triangle[0];
+ const auto& B = triangle[1];
+ const auto& C = triangle[2];
+
+ TinyMatrix<2> J;
+ for (size_t i = 0; i < 2; ++i) {
+ J(i, 0) = 0.5 * (B[i] - A[i]);
+ J(i, 1) = 0.5 * (C[i] - A[i]);
+ }
+ TinyVector s = 0.5 * (B + C);
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
- auto value = weight_list[0] * f(t(point_list[0]));
+ auto value = weight_list[0] * f(J * (point_list[0]) + s);
for (size_t i = 1; i < weight_list.size(); ++i) {
- value += weight_list[i] * f(t(point_list[i]));
+ value += weight_list[i] * f(J * (point_list[i]) + s);
}
- return t.jacobianDeterminant() * value;
+ return det(J) * value;
};
auto get_order = [&integrate, &integrate_on_rectangle](auto f, auto quadrature_formula, const double exact_value) {
diff --git a/tests/test_SimplicialGaussQuadrature.cpp b/tests/test_SimplicialGaussQuadrature.cpp
index d4e5ec7af..03b6060b5 100644
--- a/tests/test_SimplicialGaussQuadrature.cpp
+++ b/tests/test_SimplicialGaussQuadrature.cpp
@@ -526,15 +526,23 @@ TEST_CASE("SimplicialGaussQuadrature", "[analysis]")
SECTION("2D")
{
auto integrate = [](auto f, auto quadrature_formula) {
+ using R2 = TinyVector<2>;
+
+ const R2 A{-1, -1};
+ const R2 B{+1, -1};
+ const R2 C{-1, +1};
+
+ AffineTransformation<2> t{{A, B, C}};
+
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
- auto value = weight_list[0] * f(point_list[0]);
+ auto value = weight_list[0] * f(t(point_list[0]));
for (size_t i = 1; i < weight_list.size(); ++i) {
- value += weight_list[i] * f(point_list[i]);
+ value += weight_list[i] * f(t(point_list[i]));
}
- return value;
+ return t.jacobianDeterminant() * value;
};
auto integrate_on_triangle = [](auto f, auto quadrature_formula, const std::array<TinyVector<2>, 3>& triangle) {
--
GitLab