From 1fb2ac94b5b4e69a2a8e478efbf4fae43d4e9b97 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 29 Nov 2021 18:34:02 +0100
Subject: [PATCH] Allow LineTransformation in higher dimensions: ]-1,1[ -> R^d,
(d>=1)
For higher dimensions, the Jacobian determinant make now sense but the
relevant quantity (to compute integrals) is the norm of the velocity
of the curvilinear abscissa.
---
src/geometry/LineTransformation.hpp | 50 ++++++++++++++++++---
src/language/utils/IntegrateOnCells.hpp | 4 +-
tests/test_LineTransformation.cpp | 60 +++++++++++++++++++++----
3 files changed, 98 insertions(+), 16 deletions(-)
diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp
index 641821809..ae0b00a8a 100644
--- a/src/geometry/LineTransformation.hpp
+++ b/src/geometry/LineTransformation.hpp
@@ -1,16 +1,19 @@
#ifndef LINE_TRANSFORMATION_HPP
#define LINE_TRANSFORMATION_HPP
-#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-class LineTransformation
+template <size_t GivenDimension>
+class LineTransformation;
+
+template <>
+class LineTransformation<1>
{
public:
constexpr static size_t Dimension = 1;
private:
- TinyMatrix<Dimension> m_jacobian;
+ double m_jacobian;
TinyVector<Dimension> m_shift;
public:
@@ -24,14 +27,49 @@ class LineTransformation
double
jacobianDeterminant() const
{
- return m_jacobian(0, 0);
+ return m_jacobian;
+ }
+
+ PUGS_INLINE
+ LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
+ {
+ m_jacobian = 0.5 * (b[0] - a[0]);
+ m_shift = 0.5 * (a + b);
+ }
+
+ ~LineTransformation() = default;
+};
+
+template <size_t GivenDimension>
+class LineTransformation
+{
+ public:
+ constexpr static size_t Dimension = GivenDimension;
+
+ private:
+ TinyVector<Dimension> m_velocity;
+ const double m_velocity_norm;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<1>& x) const
+ {
+ return x[0] * m_velocity + m_shift;
+ }
+
+ double
+ velocityNorm() const
+ {
+ return m_velocity_norm;
}
PUGS_INLINE
LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
+ : m_velocity{0.5 * (b - a)}, m_velocity_norm{l2Norm(m_velocity)}, m_shift{0.5 * (a + b)}
{
- m_jacobian(0, 0) = 0.5 * (b[0] - a[0]);
- m_shift[0] = 0.5 * (a[0] + b[0]);
+ static_assert(Dimension > 1);
}
~LineTransformation() = default;
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 180008c15..e833fd5fb 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -140,7 +140,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
if constexpr (Dimension == 1) {
switch (cell_type[cell_id]) {
case CellType::Line: {
- const LineTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]]);
+ const LineTransformation<1> T(xr[cell_to_node[0]], xr[cell_to_node[1]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -380,7 +380,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
if constexpr (Dimension == 1) {
switch (cell_type[cell_id]) {
case CellType::Line: {
- const LineTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]]);
+ const LineTransformation<1> T(xr[cell_to_node[0]], xr[cell_to_node[1]]);
const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
diff --git a/tests/test_LineTransformation.cpp b/tests/test_LineTransformation.cpp
index f9c9620af..f7acac50f 100644
--- a/tests/test_LineTransformation.cpp
+++ b/tests/test_LineTransformation.cpp
@@ -7,16 +7,60 @@
TEST_CASE("LineTransformation", "[geometry]")
{
- using R1 = TinyVector<1>;
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
- const R1 a = 0.3;
- const R1 b = 2.7;
+ const R1 a = 0.3;
+ const R1 b = 2.7;
- const LineTransformation t(a, b);
+ const LineTransformation<1> t(a, b);
- REQUIRE(t(-1)[0] == Catch::Approx(0.3));
- REQUIRE(t(0)[0] == Catch::Approx(1.5));
- REQUIRE(t(1)[0] == Catch::Approx(2.7));
+ REQUIRE(t(-1)[0] == Catch::Approx(0.3));
+ REQUIRE(t(0)[0] == Catch::Approx(1.5));
+ REQUIRE(t(1)[0] == Catch::Approx(2.7));
- REQUIRE(t.jacobianDeterminant() == Catch::Approx(1.2));
+ REQUIRE(t.jacobianDeterminant() == Catch::Approx(1.2));
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+
+ const R2 a = {0.3, 1.2};
+ const R2 b = {2.7, 0.7};
+
+ const LineTransformation<2> t(a, b);
+
+ REQUIRE(t(-1)[0] == Catch::Approx(0.3));
+ REQUIRE(t(-1)[1] == Catch::Approx(1.2));
+ REQUIRE(t(0)[0] == Catch::Approx(1.5));
+ REQUIRE(t(0)[1] == Catch::Approx(0.95));
+ REQUIRE(t(1)[0] == Catch::Approx(2.7));
+ REQUIRE(t(1)[1] == Catch::Approx(0.7));
+
+ REQUIRE(t.velocityNorm() == Catch::Approx(l2Norm(0.5 * (b - a))));
+ }
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ const R3 a = {0.3, 1.2, -1};
+ const R3 b = {2.7, 0.7, 0.3};
+
+ const LineTransformation<3> t(a, b);
+
+ REQUIRE(t(-1)[0] == Catch::Approx(0.3));
+ REQUIRE(t(-1)[1] == Catch::Approx(1.2));
+ REQUIRE(t(-1)[2] == Catch::Approx(-1.));
+ REQUIRE(t(0)[0] == Catch::Approx(1.5));
+ REQUIRE(t(0)[1] == Catch::Approx(0.95));
+ REQUIRE(t(0)[2] == Catch::Approx(-0.35));
+ REQUIRE(t(1)[0] == Catch::Approx(2.7));
+ REQUIRE(t(1)[1] == Catch::Approx(0.7));
+ REQUIRE(t(1)[2] == Catch::Approx(0.3));
+
+ REQUIRE(t.velocityNorm() == Catch::Approx(l2Norm(0.5 * (b - a))));
+ }
}
--
GitLab