From 31b1cfd46f0bd31348872b93676940db86fa73b4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 30 Nov 2021 14:46:19 +0100
Subject: [PATCH] Allow TriangleTransformation in dimension 3
In the tests one checks that integrals are computed correctly on such surfaces.
---
src/geometry/TriangleTransformation.hpp | 48 +++++++++-
src/language/utils/IntegrateOnCells.hpp | 2 +-
tests/test_TriangleGaussQuadrature.cpp | 4 +-
tests/test_TriangleTransformation.cpp | 120 +++++++++++++++++++++---
4 files changed, 156 insertions(+), 18 deletions(-)
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index 9053912b8..9210e0bda 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -4,7 +4,11 @@
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-class TriangleTransformation
+template <size_t GivenDimension>
+class TriangleTransformation;
+
+template <>
+class TriangleTransformation<2>
{
public:
constexpr static size_t Dimension = 2;
@@ -44,4 +48,46 @@ class TriangleTransformation
~TriangleTransformation() = default;
};
+template <size_t GivenDimension>
+class TriangleTransformation
+{
+ public:
+ constexpr static size_t Dimension = GivenDimension;
+ static_assert(Dimension == 3, "Triangle transformation is defined in dimension 2 or 3");
+
+ private:
+ TinyMatrix<Dimension, 2> m_jacobian;
+ TinyVector<Dimension> m_shift;
+ double m_area_variation_norm;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<2>& x) const
+ {
+ return m_jacobian * x + m_shift;
+ }
+
+ double
+ areaVariationNorm() const
+ {
+ return m_area_variation_norm;
+ }
+
+ PUGS_INLINE
+ TriangleTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b, const TinyVector<Dimension>& c)
+ {
+ for (size_t i = 0; i < Dimension; ++i) {
+ m_jacobian(i, 0) = b[i] - a[i];
+ m_jacobian(i, 1) = c[i] - a[i];
+ }
+
+ m_shift = a;
+
+ m_area_variation_norm = l2Norm(crossProduct(b - a, c - a));
+ }
+
+ ~TriangleTransformation() = default;
+};
+
#endif // TRIANGLE_TRANSFORMATION_HPP
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index e833fd5fb..216483dc4 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -399,7 +399,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
- const TriangleTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]]);
+ const TriangleTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]]);
const auto qf = QuadratureManager::instance().getTriangleFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp
index 0bb7f34ec..3d7b5dffa 100644
--- a/tests/test_TriangleGaussQuadrature.cpp
+++ b/tests/test_TriangleGaussQuadrature.cpp
@@ -19,7 +19,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
const R2 B{+1, -1};
const R2 C{-1, +1};
- TriangleTransformation t{A, B, C};
+ TriangleTransformation<2> t{A, B, C};
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -34,7 +34,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
auto integrate_on_triangle = [](auto f, auto quadrature_formula, const TinyVector<2>& a, const TinyVector<2>& b,
const TinyVector<2>& c) {
- TriangleTransformation t(a, b, c);
+ TriangleTransformation<2> t(a, b, c);
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
diff --git a/tests/test_TriangleTransformation.cpp b/tests/test_TriangleTransformation.cpp
index a1c05aabc..99330e41d 100644
--- a/tests/test_TriangleTransformation.cpp
+++ b/tests/test_TriangleTransformation.cpp
@@ -1,31 +1,123 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/TriangleTransformation.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("TriangleTransformation", "[geometry]")
{
- using R2 = TinyVector<2>;
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
- const R2 a = {1, 2};
- const R2 b = {3, 1};
- const R2 c = {2, 5};
+ const R2 a = {1, 2};
+ const R2 b = {3, 1};
+ const R2 c = {2, 5};
- const TriangleTransformation t(a, b, c);
+ const TriangleTransformation<2> t(a, b, c);
- REQUIRE(t({0, 0})[0] == Catch::Approx(1));
- REQUIRE(t({0, 0})[1] == Catch::Approx(2));
+ REQUIRE(t({0, 0})[0] == Catch::Approx(1));
+ REQUIRE(t({0, 0})[1] == Catch::Approx(2));
- REQUIRE(t({1, 0})[0] == Catch::Approx(3));
- REQUIRE(t({1, 0})[1] == Catch::Approx(1));
+ REQUIRE(t({1, 0})[0] == Catch::Approx(3));
+ REQUIRE(t({1, 0})[1] == Catch::Approx(1));
- REQUIRE(t({0, 1})[0] == Catch::Approx(2));
- REQUIRE(t({0, 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));
+ REQUIRE(t.jacobianDeterminant() == Catch::Approx(7));
+ }
+
+ SECTION("3D")
+ {
+ SECTION("Data")
+ {
+ using R3 = TinyVector<3>;
+
+ const R3 a = {1, 2, 2};
+ const R3 b = {4, 1, 3};
+ const R3 c = {2, 5, 1};
+
+ const TriangleTransformation<3> t(a, b, c);
+
+ REQUIRE(t({0, 0})[0] == Catch::Approx(1));
+ REQUIRE(t({0, 0})[1] == Catch::Approx(2));
+ REQUIRE(t({0, 0})[2] == Catch::Approx(2));
+
+ REQUIRE(t({1, 0})[0] == Catch::Approx(4));
+ REQUIRE(t({1, 0})[1] == Catch::Approx(1));
+ REQUIRE(t({1, 0})[2] == Catch::Approx(3));
+
+ REQUIRE(t({0, 1})[0] == Catch::Approx(2));
+ REQUIRE(t({0, 1})[1] == Catch::Approx(5));
+ REQUIRE(t({0, 1})[2] == Catch::Approx(1));
+
+ REQUIRE(t({1. / 3, 1. / 3})[0] == Catch::Approx(7. / 3));
+ REQUIRE(t({1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
+ REQUIRE(t({1. / 3, 1. / 3})[2] == Catch::Approx(2));
+
+ REQUIRE(t.areaVariationNorm() == Catch::Approx(2 * std::sqrt(30)));
+ }
+
+ SECTION("Area modulus")
+ {
+ using R3 = TinyVector<3>;
+
+ const R3 a_hat = {0, 0, 0};
+ const R3 b_hat = {1, 0, 0};
+ const R3 c_hat = {0, 1, 0};
+
+ const double theta = 2.3;
+ TinyMatrix<3> r_theta(1, 0, 0, //
+ 0, std::cos(theta), std::sin(theta), //
+ 0, -std::sin(theta), std::cos(theta));
+
+ const double phi = 0.7;
+ TinyMatrix<3> r_phi(std::cos(phi), std::sin(phi), 0, //
+ -std::sin(phi), std::cos(phi), 0, //
+ 0, 0, 1);
+
+ const R3 a = r_phi * r_theta * a_hat;
+ const R3 b = r_phi * r_theta * b_hat;
+ const R3 c = r_phi * r_theta * c_hat;
+
+ const TriangleTransformation<3> t(a, b, c);
+
+ // The triangle (a,b,c) is just a rotation of the initial one
+ REQUIRE(t.areaVariationNorm() == Catch::Approx(1));
+ }
+
+ SECTION("Polynomial integral")
+ {
+ using R3 = TinyVector<3>;
+
+ const R3 a = {1, 2, 2};
+ const R3 b = {4, 1, 3};
+ const R3 c = {2, 5, 1};
+
+ const TriangleTransformation<3> t(a, b, c);
+
+ 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 - 3 * y * y + y + 2 * z * z - 0.5 * z + 2;
+ };
+
+ 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.areaVariationNorm() * p(t(qf.point(i)));
+ }
+
+ REQUIRE(sum == Catch::Approx(39.2534499545369));
+ }
+ }
}
--
GitLab