diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp index 641821809f583e007851c259e959584fbcc2e48c..ae0b00a8a1dce7e4b6b6bb13a737ee1114a202f1 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 180008c155218514949038f4619b21ef7192cd8f..e833fd5fb1fd04990fb169687776be7ee0d311cb 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 f9c9620af26851da7fd3724a6d4f8c39436d301f..f7acac50fdb409bc2738d156fc1b9833e8ffff0e 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)))); + } }