diff --git a/src/geometry/AffineTransformation.hpp b/src/geometry/AffineTransformation.hpp
deleted file mode 100644
index ed5a5509b39a6439eea197058d1d69df439b6768..0000000000000000000000000000000000000000
--- a/src/geometry/AffineTransformation.hpp
+++ /dev/null
@@ -1,78 +0,0 @@
-#ifndef AFFINE_TRANSFORMATION_HPP
-#define AFFINE_TRANSFORMATION_HPP
-
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-
-#include <array>
-
-template <size_t Dimension>
-class AffineTransformation
-{
- private:
- constexpr static size_t NumberOfPoints = Dimension + 1;
-
- TinyMatrix<Dimension> m_jacobian;
- TinyVector<Dimension> m_shift;
- double m_jacobian_determinant;
-
- public:
- PUGS_INLINE
- TinyVector<Dimension>
- operator()(const TinyVector<Dimension>& x) const
- {
- return m_jacobian * x + m_shift;
- }
-
- double
- jacobianDeterminant() const
- {
- return m_jacobian_determinant;
- }
-
- PUGS_INLINE
- AffineTransformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
- {
- static_assert(Dimension >= 1 and Dimension <= 3, "invalid dimension");
- if constexpr (Dimension == 1) {
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
-
- m_jacobian(0, 0) = 0.5 * (b[0] - a[0]);
-
- m_shift[0] = 0.5 * (a[0] + b[0]);
- } else if constexpr (Dimension == 2) {
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
-
- 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;
- } else {
- static_assert(Dimension == 3);
-
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
- const TinyVector<Dimension>& d = image_nodes[3];
-
- 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_jacobian(i, 2) = d[i] - a[i];
- }
-
- m_shift = a;
- }
-
- m_jacobian_determinant = det(m_jacobian);
- }
-
- ~AffineTransformation() = default;
-};
-
-#endif // AFFINE_TRANSFORMATION_HPP
diff --git a/src/geometry/CubeTransformation.hpp b/src/geometry/CubeTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7bc12c324e4a55750b41b1584633cc789f8a7b2d
--- /dev/null
+++ b/src/geometry/CubeTransformation.hpp
@@ -0,0 +1,79 @@
+#ifndef CUBE_TRANSFORMATION_HPP
+#define CUBE_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+class CubeTransformation
+{
+ public:
+ constexpr static size_t Dimension = 3;
+ constexpr static size_t NumberOfPoints = 8;
+
+ private:
+ TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<Dimension>& x) const
+ {
+ const TinyVector<NumberOfPoints - 1> X =
+ {x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2], x[0] * x[1] * x[2]};
+ return m_coefficients * X + m_shift;
+ }
+
+ PUGS_INLINE double
+ jacobianDeterminant(const TinyVector<Dimension>& X) const
+ {
+ static_assert(Dimension == 3, "invalid dimension");
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+ const double& z = X[2];
+
+ const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * y + T(0, 5) * z + T(0, 6) * y * z,
+ T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * z,
+ T(0, 2) + T(0, 4) * y + T(0, 5) * x + T(0, 6) * x * y,
+ //
+ T(1, 0) + T(1, 3) * y + T(1, 5) * z + T(1, 6) * y * z,
+ T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * z,
+ T(1, 2) + T(1, 4) * y + T(1, 5) * x + T(1, 6) * x * y,
+ //
+ T(2, 0) + T(2, 3) * y + T(2, 5) * z + T(2, 6) * y * z,
+ T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * z,
+ T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
+
+ return det(J);
+ }
+
+ PUGS_INLINE
+ CubeTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d,
+ const TinyVector<Dimension>& e,
+ const TinyVector<Dimension>& f,
+ const TinyVector<Dimension>& g,
+ const TinyVector<Dimension>& h)
+ {
+ for (size_t i = 0; i < Dimension; ++i) {
+ m_coefficients(i, 0) = 0.125 * (-a[i] + b[i] + c[i] - d[i] - e[i] + f[i] + g[i] - h[i]);
+ m_coefficients(i, 1) = 0.125 * (-a[i] - b[i] + c[i] + d[i] - e[i] - f[i] + g[i] + h[i]);
+ m_coefficients(i, 2) = 0.125 * (-a[i] - b[i] - c[i] - d[i] + e[i] + f[i] + g[i] + h[i]);
+ m_coefficients(i, 3) = 0.125 * (+a[i] - b[i] + c[i] - d[i] + e[i] - f[i] + g[i] - h[i]);
+ m_coefficients(i, 4) = 0.125 * (+a[i] + b[i] - c[i] - d[i] - e[i] - f[i] + g[i] + h[i]);
+ m_coefficients(i, 5) = 0.125 * (+a[i] - b[i] - c[i] + d[i] - e[i] + f[i] + g[i] - h[i]);
+ m_coefficients(i, 6) = 0.125 * (-a[i] + b[i] - c[i] + d[i] + e[i] - f[i] + g[i] - h[i]);
+
+ m_shift[i] = 0.125 * (a[i] + b[i] + c[i] + d[i] + e[i] + f[i] + g[i] + h[i]);
+ }
+ }
+
+ ~CubeTransformation() = default;
+};
+
+#endif // CUBE_TRANSFORMATION_HPP
diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..641821809f583e007851c259e959584fbcc2e48c
--- /dev/null
+++ b/src/geometry/LineTransformation.hpp
@@ -0,0 +1,40 @@
+#ifndef LINE_TRANSFORMATION_HPP
+#define LINE_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+class LineTransformation
+{
+ public:
+ constexpr static size_t Dimension = 1;
+
+ private:
+ TinyMatrix<Dimension> m_jacobian;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<1>& x) const
+ {
+ return m_jacobian * x + m_shift;
+ }
+
+ double
+ jacobianDeterminant() const
+ {
+ return m_jacobian(0, 0);
+ }
+
+ PUGS_INLINE
+ LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
+ {
+ m_jacobian(0, 0) = 0.5 * (b[0] - a[0]);
+ m_shift[0] = 0.5 * (a[0] + b[0]);
+ }
+
+ ~LineTransformation() = default;
+};
+
+#endif // LINE_TRANSFORMATION_HPP
diff --git a/src/geometry/PrismTransformation.hpp b/src/geometry/PrismTransformation.hpp
index 730710bba373584d4c21c8e3dfbedf20a62e96b9..5e6e509ed0c92023b81ea13a3e686392dec52fcd 100644
--- a/src/geometry/PrismTransformation.hpp
+++ b/src/geometry/PrismTransformation.hpp
@@ -4,8 +4,6 @@
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-#include <array>
-
class PrismTransformation
{
private:
@@ -47,17 +45,15 @@ class PrismTransformation
}
PUGS_INLINE
- PrismTransformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
+ PrismTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d,
+ const TinyVector<Dimension>& e,
+ const TinyVector<Dimension>& f)
{
static_assert(Dimension == 3, "invalid dimension");
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
- const TinyVector<Dimension>& d = image_nodes[3];
- const TinyVector<Dimension>& e = image_nodes[4];
- const TinyVector<Dimension>& f = image_nodes[5];
-
for (size_t i = 0; i < Dimension; ++i) {
m_coefficients(i, 0) = 0.5 * (b[i] - a[i] + e[i] - d[i]);
m_coefficients(i, 1) = 0.5 * (c[i] + f[i] - a[i] - d[i]);
diff --git a/src/geometry/PyramidTransformation.hpp b/src/geometry/PyramidTransformation.hpp
index 19e7822ba3bc3e01da7be2840522e24707b68007..2fa3eaa8d8b010dda5ac3d591b1183db707ea3f8 100644
--- a/src/geometry/PyramidTransformation.hpp
+++ b/src/geometry/PyramidTransformation.hpp
@@ -4,8 +4,6 @@
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-#include <array>
-
class PyramidTransformation
{
private:
@@ -46,16 +44,14 @@ class PyramidTransformation
}
PUGS_INLINE
- PyramidTransformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
+ PyramidTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d,
+ const TinyVector<Dimension>& e)
{
static_assert(Dimension == 3, "invalid dimension");
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
- const TinyVector<Dimension>& d = image_nodes[3];
- const TinyVector<Dimension>& e = image_nodes[4];
-
m_shift = 0.25 * (a + b + c + d);
for (size_t i = 0; i < Dimension; ++i) {
diff --git a/src/geometry/Q1Transformation.hpp b/src/geometry/Q1Transformation.hpp
deleted file mode 100644
index 55e87bcbb11fede5993e8f0bb8f9a19881b76352..0000000000000000000000000000000000000000
--- a/src/geometry/Q1Transformation.hpp
+++ /dev/null
@@ -1,128 +0,0 @@
-#ifndef Q1_TRANSFORMATION_HPP
-#define Q1_TRANSFORMATION_HPP
-
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-
-#include <array>
-
-template <size_t Dimension>
-class Q1Transformation
-{
- private:
- constexpr static size_t NumberOfPoints = 1 << Dimension;
-
- TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
- TinyVector<Dimension> m_shift;
-
- public:
- PUGS_INLINE
- TinyVector<Dimension>
- operator()(const TinyVector<Dimension>& x) const
- {
- if constexpr (Dimension == 1) {
- return m_coefficients * x + m_shift;
- } else if constexpr (Dimension == 2) {
- const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
- return m_coefficients * X + m_shift;
- } else {
- static_assert(Dimension == 3, "invalid dimension");
- const TinyVector<NumberOfPoints - 1> X =
- {x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2], x[0] * x[1] * x[2]};
- return m_coefficients * X + m_shift;
- }
- }
-
- PUGS_INLINE double
- jacobianDeterminant([[maybe_unused]] const TinyVector<Dimension>& X) const
- {
- if constexpr (Dimension == 1) {
- return m_coefficients(0, 0);
- } else if constexpr (Dimension == 2) {
- const auto& T = m_coefficients;
- const double& x = X[0];
- const double& y = X[1];
-
- const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, //
- T(0, 1) + T(0, 2) * x,
- //
- T(1, 0) + T(1, 2) * y, //
- T(1, 1) + T(1, 2) * x};
- return det(J);
- } else {
- static_assert(Dimension == 3, "invalid dimension");
- const auto& T = m_coefficients;
- const double& x = X[0];
- const double& y = X[1];
- const double& z = X[2];
-
- const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * y + T(0, 5) * z + T(0, 6) * y * z,
- T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * z,
- T(0, 2) + T(0, 4) * y + T(0, 5) * x + T(0, 6) * x * y,
- //
- T(1, 0) + T(1, 3) * y + T(1, 5) * z + T(1, 6) * y * z,
- T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * z,
- T(1, 2) + T(1, 4) * y + T(1, 5) * x + T(1, 6) * x * y,
- //
- T(2, 0) + T(2, 3) * y + T(2, 5) * z + T(2, 6) * y * z,
- T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * z,
- T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
-
- return det(J);
- }
- }
-
- PUGS_INLINE
- Q1Transformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
- {
- static_assert(Dimension >= 1 and Dimension <= 3, "invalid dimension");
- if constexpr (Dimension == 1) {
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
-
- m_coefficients(0, 0) = 0.5 * (b[0] - a[0]);
-
- m_shift[0] = 0.5 * (a[0] + b[0]);
- } else if constexpr (Dimension == 2) {
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
- const TinyVector<Dimension>& d = image_nodes[3];
-
- for (size_t i = 0; i < Dimension; ++i) {
- m_coefficients(i, 0) = 0.25 * (-a[i] + b[i] + c[i] - d[i]);
- m_coefficients(i, 1) = 0.25 * (-a[i] - b[i] + c[i] + d[i]);
- m_coefficients(i, 2) = 0.25 * (+a[i] - b[i] + c[i] - d[i]);
-
- m_shift[i] = 0.25 * (a[i] + b[i] + c[i] + d[i]);
- }
- } else {
- static_assert(Dimension == 3);
-
- const TinyVector<Dimension>& a = image_nodes[0];
- const TinyVector<Dimension>& b = image_nodes[1];
- const TinyVector<Dimension>& c = image_nodes[2];
- const TinyVector<Dimension>& d = image_nodes[3];
- const TinyVector<Dimension>& e = image_nodes[4];
- const TinyVector<Dimension>& f = image_nodes[5];
- const TinyVector<Dimension>& g = image_nodes[6];
- const TinyVector<Dimension>& h = image_nodes[7];
-
- for (size_t i = 0; i < Dimension; ++i) {
- m_coefficients(i, 0) = 0.125 * (-a[i] + b[i] + c[i] - d[i] - e[i] + f[i] + g[i] - h[i]);
- m_coefficients(i, 1) = 0.125 * (-a[i] - b[i] + c[i] + d[i] - e[i] - f[i] + g[i] + h[i]);
- m_coefficients(i, 2) = 0.125 * (-a[i] - b[i] - c[i] - d[i] + e[i] + f[i] + g[i] + h[i]);
- m_coefficients(i, 3) = 0.125 * (+a[i] - b[i] + c[i] - d[i] + e[i] - f[i] + g[i] - h[i]);
- m_coefficients(i, 4) = 0.125 * (+a[i] + b[i] - c[i] - d[i] - e[i] - f[i] + g[i] + h[i]);
- m_coefficients(i, 5) = 0.125 * (+a[i] - b[i] - c[i] + d[i] - e[i] + f[i] + g[i] - h[i]);
- m_coefficients(i, 6) = 0.125 * (-a[i] + b[i] - c[i] + d[i] + e[i] - f[i] + g[i] - h[i]);
-
- m_shift[i] = 0.125 * (a[i] + b[i] + c[i] + d[i] + e[i] + f[i] + g[i] + h[i]);
- }
- }
- }
-
- ~Q1Transformation() = default;
-};
-
-#endif // Q1_TRANSFORMATION_HPP
diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d2906bb514fe90137f7270404c2f5b5f54d24a2
--- /dev/null
+++ b/src/geometry/SquareTransformation.hpp
@@ -0,0 +1,61 @@
+#ifndef SQUARE_TRANSFORMATION_HPP
+#define SQUARE_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+class SquareTransformation
+{
+ public:
+ constexpr static size_t Dimension = 2;
+ constexpr static size_t NumberOfPoints = 4;
+
+ private:
+ TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<2>& x) const
+ {
+ const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+ return m_coefficients * X + m_shift;
+ }
+
+ PUGS_INLINE double
+ jacobianDeterminant(const TinyVector<Dimension>& X) const
+ {
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+
+ const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, //
+ T(0, 1) + T(0, 2) * x,
+ //
+ T(1, 0) + T(1, 2) * y, //
+ T(1, 1) + T(1, 2) * x};
+ return det(J);
+ }
+
+ PUGS_INLINE
+ SquareTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d)
+ {
+ for (size_t i = 0; i < Dimension; ++i) {
+ m_coefficients(i, 0) = 0.25 * (-a[i] + b[i] + c[i] - d[i]);
+ m_coefficients(i, 1) = 0.25 * (-a[i] - b[i] + c[i] + d[i]);
+ m_coefficients(i, 2) = 0.25 * (+a[i] - b[i] + c[i] - d[i]);
+
+ m_shift[i] = 0.25 * (a[i] + b[i] + c[i] + d[i]);
+ }
+ }
+
+ ~SquareTransformation() = default;
+};
+
+#endif // SQUARE_TRANSFORMATION_HPP
diff --git a/src/geometry/TetrahedronTransformation.hpp b/src/geometry/TetrahedronTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..64244266957a350e8be8f78c2eab3d5fead629e1
--- /dev/null
+++ b/src/geometry/TetrahedronTransformation.hpp
@@ -0,0 +1,53 @@
+#ifndef TETRAHEDRON_TRANSFORMATION_HPP
+#define TETRAHEDRON_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+class TetrahedronTransformation
+{
+ public:
+ constexpr static size_t Dimension = 3;
+
+ private:
+ constexpr static size_t NumberOfPoints = Dimension + 1;
+
+ TinyMatrix<Dimension> m_jacobian;
+ TinyVector<Dimension> m_shift;
+ double m_jacobian_determinant;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<Dimension>& x) const
+ {
+ return m_jacobian * x + m_shift;
+ }
+
+ double
+ jacobianDeterminant() const
+ {
+ return m_jacobian_determinant;
+ }
+
+ PUGS_INLINE
+ TetrahedronTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d)
+ {
+ 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_jacobian(i, 2) = d[i] - a[i];
+ }
+
+ m_shift = a;
+
+ m_jacobian_determinant = det(m_jacobian);
+ }
+
+ ~TetrahedronTransformation() = default;
+};
+
+#endif // TETRAHEDRON_TRANSFORMATION_HPP
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9053912b8057c6d1c2e0ba8227a9f889bb1fc8f4
--- /dev/null
+++ b/src/geometry/TriangleTransformation.hpp
@@ -0,0 +1,47 @@
+#ifndef TRIANGLE_TRANSFORMATION_HPP
+#define TRIANGLE_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+class TriangleTransformation
+{
+ public:
+ constexpr static size_t Dimension = 2;
+
+ private:
+ TinyMatrix<Dimension> m_jacobian;
+ TinyVector<Dimension> m_shift;
+ double m_jacobian_determinant;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<2>& x) const
+ {
+ return m_jacobian * x + m_shift;
+ }
+
+ double
+ jacobianDeterminant() const
+ {
+ return m_jacobian_determinant;
+ }
+
+ 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_jacobian_determinant = det(m_jacobian);
+ }
+
+ ~TriangleTransformation() = default;
+};
+
+#endif // TRIANGLE_TRANSFORMATION_HPP
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index d4463da9647efc1877ba6df38e4042d17a76f168..7145c662bb5e294f4f1235bc51665f020fd8bd04 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -3,10 +3,13 @@
#include <analysis/IQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
-#include <geometry/AffineTransformation.hpp>
+#include <geometry/CubeTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
#include <geometry/PrismTransformation.hpp>
#include <geometry/PyramidTransformation.hpp>
-#include <geometry/Q1Transformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
#include <language/utils/PugsFunctionAdapter.hpp>
#include <mesh/CellType.hpp>
#include <mesh/Connectivity.hpp>
@@ -137,7 +140,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
if constexpr (Dimension == 1) {
switch (cell_type[cell_id]) {
case CellType::Line: {
- const AffineTransformation<1> T({xr[cell_to_node[0]], xr[cell_to_node[1]]});
+ const LineTransformation 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);
@@ -155,8 +158,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
- const Q1Transformation<2> T(
- {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[2]]});
+ const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[2]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -167,8 +170,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Quadrangle: {
- const Q1Transformation<2> T(
- {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]]});
+ const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -188,10 +191,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
switch (cell_type[cell_id]) {
case CellType::Tetrahedron: {
- const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[2]], //
- xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[3]],
- xr[cell_to_node[3]]});
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[3]],
+ xr[cell_to_node[3]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -203,10 +206,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
case CellType::Pyramid: {
if (cell_to_node.size() == 5) {
- const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
- xr[cell_to_node[3]], //
- xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
- xr[cell_to_node[4]]});
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -223,10 +226,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
case CellType::Diamond: {
if (cell_to_node.size() == 5) {
{ // top tetrahedron
- const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
- xr[cell_to_node[3]], //
- xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
- xr[cell_to_node[4]]});
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[3]], //
+ xr[cell_to_node[4]], xr[cell_to_node[4]], xr[cell_to_node[4]],
+ xr[cell_to_node[4]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -236,10 +239,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
{ // bottom tetrahedron
- const Q1Transformation<3> T1({xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[2]],
- xr[cell_to_node[1]], //
- xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
- xr[cell_to_node[0]]});
+ const CubeTransformation T1(xr[cell_to_node[3]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -250,10 +253,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
} else if (cell_to_node.size() == 6) {
{ // top pyramid
- const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
- xr[cell_to_node[4]], //
- xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
- xr[cell_to_node[5]]});
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -263,10 +266,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
{ // bottom pyramid
- const Q1Transformation<3> T1({xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
- xr[cell_to_node[1]], //
- xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
- xr[cell_to_node[0]]});
+ const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -282,10 +285,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Prism: {
- const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], //
- xr[cell_to_node[2]], xr[cell_to_node[2]], //
- xr[cell_to_node[3]], xr[cell_to_node[4]], //
- xr[cell_to_node[5]], xr[cell_to_node[5]]});
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -296,9 +299,9 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Hexahedron: {
- const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]],
- xr[cell_to_node[6]], xr[cell_to_node[7]]});
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]], xr[cell_to_node[6]],
+ xr[cell_to_node[7]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -377,7 +380,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
if constexpr (Dimension == 1) {
switch (cell_type[cell_id]) {
case CellType::Line: {
- const AffineTransformation<1> T({xr[cell_to_node[0]], xr[cell_to_node[1]]});
+ const LineTransformation 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) {
@@ -396,7 +399,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
- const AffineTransformation<2> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]]});
+ const TriangleTransformation 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) {
@@ -408,8 +411,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Quadrangle: {
- const Q1Transformation<2> T(
- {xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]]});
+ const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]]);
const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -430,8 +433,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
switch (cell_type[cell_id]) {
case CellType::Tetrahedron: {
- const AffineTransformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], //
- xr[cell_to_node[2]], xr[cell_to_node[3]]});
+ const TetrahedronTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], //
+ xr[cell_to_node[2]], xr[cell_to_node[3]]);
const auto qf = QuadratureManager::instance().getTetrahedronFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -444,8 +447,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
case CellType::Pyramid: {
if (cell_to_node.size() == 5) {
- const PyramidTransformation T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
- xr[cell_to_node[3]], xr[cell_to_node[4]]});
+ const PyramidTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], //
+ xr[cell_to_node[3]], xr[cell_to_node[4]]);
const auto qf = QuadratureManager::instance().getPyramidFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -463,8 +466,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
if (cell_to_node.size() == 5) {
const auto qf = QuadratureManager::instance().getTetrahedronFormula(quadrature_descriptor);
{ // top tetrahedron
- const AffineTransformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]], //
- xr[cell_to_node[4]]});
+ const TetrahedronTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]], //
+ xr[cell_to_node[4]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -474,8 +477,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
{ // bottom tetrahedron
- const AffineTransformation<3> T1({xr[cell_to_node[3]], xr[cell_to_node[2]], xr[cell_to_node[1]], //
- xr[cell_to_node[0]]});
+ const TetrahedronTransformation T1(xr[cell_to_node[3]], xr[cell_to_node[2]], xr[cell_to_node[1]], //
+ xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -487,10 +490,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
} else if (cell_to_node.size() == 6) {
const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
{ // top pyramid
- const Q1Transformation<3> T0({xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
- xr[cell_to_node[4]], //
- xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
- xr[cell_to_node[5]]});
+ const CubeTransformation T0(xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], //
+ xr[cell_to_node[5]], xr[cell_to_node[5]], xr[cell_to_node[5]],
+ xr[cell_to_node[5]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -500,10 +503,10 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
}
}
{ // bottom pyramid
- const Q1Transformation<3> T1({xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
- xr[cell_to_node[1]], //
- xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
- xr[cell_to_node[0]]});
+ const CubeTransformation T1(xr[cell_to_node[4]], xr[cell_to_node[3]], xr[cell_to_node[2]],
+ xr[cell_to_node[1]], //
+ xr[cell_to_node[0]], xr[cell_to_node[0]], xr[cell_to_node[0]],
+ xr[cell_to_node[0]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -519,8 +522,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Prism: {
- const PrismTransformation T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]]});
+ const PrismTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]]);
const auto qf = QuadratureManager::instance().getPrismFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -532,9 +535,9 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Hexahedron: {
- const Q1Transformation<3> T({xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[3]], xr[cell_to_node[4]], xr[cell_to_node[5]],
- xr[cell_to_node[6]], xr[cell_to_node[7]]});
+ const CubeTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]], xr[cell_to_node[3]],
+ xr[cell_to_node[4]], xr[cell_to_node[5]], xr[cell_to_node[6]],
+ xr[cell_to_node[7]]);
const auto qf = QuadratureManager::instance().getCubeFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index edd77f19444d13ec9024241d686eaa66085f6ada..0994985c3e5fa5e30198e1374aee1c0f8713025c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -9,7 +9,6 @@ add_executable (unit_tests
test_AffectationProcessor.cpp
test_AffectationToStringProcessor.cpp
test_AffectationToTupleProcessor.cpp
- test_AffineTransformation.cpp
test_Array.cpp
test_ArraySubscriptProcessor.cpp
test_ASTBuilder.cpp
@@ -62,6 +61,7 @@ add_executable (unit_tests
test_CRSMatrix.cpp
test_CRSMatrixDescriptor.cpp
test_CubeGaussQuadrature.cpp
+ test_CubeTransformation.cpp
test_DataVariant.cpp
test_Demangle.cpp
test_DiscreteFunctionDescriptorP0.cpp
@@ -87,6 +87,7 @@ add_executable (unit_tests
test_ItemType.cpp
test_LinearSolver.cpp
test_LinearSolverOptions.cpp
+ test_LineTransformation.cpp
test_ListAffectationProcessor.cpp
test_MathModule.cpp
test_NameProcessor.cpp
@@ -101,20 +102,22 @@ add_executable (unit_tests
test_PugsUtils.cpp
test_PyramidGaussQuadrature.cpp
test_PyramidTransformation.cpp
- test_Q1Transformation.cpp
test_RevisionInfo.cpp
test_SmallArray.cpp
test_SmallVector.cpp
test_SquareGaussQuadrature.cpp
+ test_SquareTransformation.cpp
test_SymbolTable.cpp
test_Table.cpp
test_TetrahedronGaussQuadrature.cpp
+ test_TetrahedronTransformation.cpp
test_TensorialGaussLegendreQuadrature.cpp
test_TensorialGaussLobattoQuadrature.cpp
test_Timer.cpp
test_TinyMatrix.cpp
test_TinyVector.cpp
test_TriangleGaussQuadrature.cpp
+ test_TriangleTransformation.cpp
test_TupleToVectorProcessor.cpp
test_UnaryExpressionProcessor.cpp
test_UnaryOperatorMangler.cpp
diff --git a/tests/test_AffineTransformation.cpp b/tests/test_AffineTransformation.cpp
deleted file mode 100644
index ec1dfa482434694f39c895ec3e9bf2e00fbba28e..0000000000000000000000000000000000000000
--- a/tests/test_AffineTransformation.cpp
+++ /dev/null
@@ -1,88 +0,0 @@
-#include <catch2/catch_approx.hpp>
-#include <catch2/catch_test_macros.hpp>
-
-#include <geometry/AffineTransformation.hpp>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("AffineTransformation", "[geometry]")
-{
- SECTION("1D")
- {
- using R1 = TinyVector<1>;
-
- const R1 a = 0.3;
- const R1 b = 2.7;
-
- const std::array points = {a, b};
-
- const AffineTransformation<1> t(points);
-
- 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));
- }
-
- SECTION("2D")
- {
- using R2 = TinyVector<2>;
-
- const R2 a = {1, 2};
- const R2 b = {3, 1};
- const R2 c = {2, 5};
-
- const std::array points = {a, b, c};
- const AffineTransformation<2> t(points);
-
- 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({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.jacobianDeterminant() == Catch::Approx(7));
- }
-
- SECTION("3D")
- {
- using R3 = TinyVector<3>;
-
- const R3 a = {1, 2, 1};
- const R3 b = {3, 1, 3};
- const R3 c = {2, 5, 2};
- const R3 d = {2, 3, 4};
-
- const std::array points = {a, b, c, d};
- const AffineTransformation<3> t(points);
-
- REQUIRE(t({0, 0, 0})[0] == Catch::Approx(1));
- REQUIRE(t({0, 0, 0})[1] == Catch::Approx(2));
- REQUIRE(t({0, 0, 0})[2] == Catch::Approx(1));
-
- REQUIRE(t({1, 0, 0})[0] == Catch::Approx(3));
- REQUIRE(t({1, 0, 0})[1] == Catch::Approx(1));
- REQUIRE(t({1, 0, 0})[2] == Catch::Approx(3));
-
- REQUIRE(t({0, 1, 0})[0] == Catch::Approx(2));
- REQUIRE(t({0, 1, 0})[1] == Catch::Approx(5));
- REQUIRE(t({0, 1, 0})[2] == Catch::Approx(2));
-
- REQUIRE(t({0, 0, 1})[0] == Catch::Approx(2));
- REQUIRE(t({0, 0, 1})[1] == Catch::Approx(3));
- REQUIRE(t({0, 0, 1})[2] == Catch::Approx(4));
-
- REQUIRE(t({0.25, 0.25, 0.25})[0] == Catch::Approx(2));
- REQUIRE(t({0.25, 0.25, 0.25})[1] == Catch::Approx(11. / 4));
- REQUIRE(t({0.25, 0.25, 0.25})[2] == Catch::Approx(2.5));
-
- REQUIRE(t.jacobianDeterminant() == Catch::Approx(14));
- }
-}
diff --git a/tests/test_CubeTransformation.cpp b/tests/test_CubeTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2de108f9c96c5c3659c20b4c3e7d314d5342280c
--- /dev/null
+++ b/tests/test_CubeTransformation.cpp
@@ -0,0 +1,117 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("CubeTransformation", "[geometry]")
+{
+ using R3 = TinyVector<3>;
+
+ const R3 a_hat = {-1, -1, -1};
+ const R3 b_hat = {+1, -1, -1};
+ const R3 c_hat = {+1, +1, -1};
+ const R3 d_hat = {-1, +1, -1};
+ const R3 e_hat = {-1, -1, +1};
+ const R3 f_hat = {+1, -1, +1};
+ const R3 g_hat = {+1, +1, +1};
+ const R3 h_hat = {-1, +1, +1};
+
+ const R3 m_hat = zero;
+
+ const R3 a = {0, 0, 0};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {0, 3, 1};
+ const R3 e = {1, 2, 5};
+ const R3 f = {3, 1, 7};
+ const R3 g = {2, 5, 5};
+ const R3 h = {0, 3, 6};
+
+ const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
+
+ const CubeTransformation t(a, b, c, d, e, f, g, h);
+
+ SECTION("values")
+ {
+ REQUIRE(l2Norm(t(a_hat) - a) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(b_hat) - b) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(c_hat) - c) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(d_hat) - d) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(e_hat) - e) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(g_hat) - g) == Catch::Approx(0));
+ REQUIRE(l2Norm(t(h_hat) - h) == Catch::Approx(0));
+
+ REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0));
+ }
+
+ SECTION("Jacobian determinant")
+ {
+ SECTION("at points")
+ {
+ auto detJ = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+
+ return -(3 * x * y * y * z + 9 * y * y * z - 6 * x * x * y * z + 6 * x * y * z - 96 * y * z + 14 * x * x * z -
+ 49 * x * z + 119 * z - 21 * x * y * y + 21 * y * y + 90 * x * y - 10 * y + 40 * x * x - 109 * x -
+ 491) /
+ 128;
+ };
+
+ REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
+ REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
+ REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
+ REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
+ REQUIRE(t.jacobianDeterminant(e_hat) == Catch::Approx(detJ(e_hat)));
+ REQUIRE(t.jacobianDeterminant(f_hat) == Catch::Approx(detJ(f_hat)));
+ REQUIRE(t.jacobianDeterminant(g_hat) == Catch::Approx(detJ(g_hat)));
+ REQUIRE(t.jacobianDeterminant(h_hat) == Catch::Approx(detJ(h_hat)));
+
+ REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+ }
+
+ SECTION("Gauss order 3")
+ {
+ // The jacobian determinant is of maximal degree 2 in each variable
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(2));
+
+ double volume = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ }
+
+ // 353. / 12 is actually the volume of the hexahedron
+ REQUIRE(volume == Catch::Approx(353. / 12));
+ }
+
+ SECTION("Gauss order 4")
+ {
+ 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 * y * z + z * z + 2 * x - 3 * y + 0.5 * z + 3;
+ };
+
+ // Jacbian determinant is a degree 2 polynomial, so the
+ // following formula is required to reach exactness
+ const QuadratureFormula<3>& gauss =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ }
+
+ REQUIRE(integral == Catch::Approx(8795. / 72));
+ }
+ }
+}
diff --git a/tests/test_LineTransformation.cpp b/tests/test_LineTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9c9620af26851da7fd3724a6d4f8c39436d301f
--- /dev/null
+++ b/tests/test_LineTransformation.cpp
@@ -0,0 +1,22 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <geometry/LineTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("LineTransformation", "[geometry]")
+{
+ using R1 = TinyVector<1>;
+
+ const R1 a = 0.3;
+ const R1 b = 2.7;
+
+ const LineTransformation 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.jacobianDeterminant() == Catch::Approx(1.2));
+}
diff --git a/tests/test_PrismGaussQuadrature.cpp b/tests/test_PrismGaussQuadrature.cpp
index e1956601c591922fcc4e57acedb12b7bc351f903..ceff287e5ac5c823c9f427ba6801babfeaa47481 100644
--- a/tests/test_PrismGaussQuadrature.cpp
+++ b/tests/test_PrismGaussQuadrature.cpp
@@ -5,7 +5,7 @@
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/PrismGaussQuadrature.hpp>
#include <analysis/QuadratureManager.hpp>
-#include <geometry/AffineTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
diff --git a/tests/test_PrismTransformation.cpp b/tests/test_PrismTransformation.cpp
index 5b19f9b08ca5631370eed4a4787e9c521a8f7ccf..e170da3c05cf2643ab22958509b2cf62039f3fa1 100644
--- a/tests/test_PrismTransformation.cpp
+++ b/tests/test_PrismTransformation.cpp
@@ -27,8 +27,7 @@ TEST_CASE("PrismTransformation", "[geometry]")
const R3 e = {1, 2, 5};
const R3 f = {3, 1, 7};
- const std::array points = {a, b, c, d, e, f};
- const PrismTransformation t(points);
+ const PrismTransformation t(a, b, c, d, e, f);
SECTION("points")
{
diff --git a/tests/test_PyramidGaussQuadrature.cpp b/tests/test_PyramidGaussQuadrature.cpp
index 9c8bfc9cc301d64d486fa617b7630ff88f012df2..e648c96fb13df144945e7c73cdf381545ff1169b 100644
--- a/tests/test_PyramidGaussQuadrature.cpp
+++ b/tests/test_PyramidGaussQuadrature.cpp
@@ -5,7 +5,7 @@
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/PyramidGaussQuadrature.hpp>
#include <analysis/QuadratureManager.hpp>
-#include <geometry/AffineTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
diff --git a/tests/test_PyramidTransformation.cpp b/tests/test_PyramidTransformation.cpp
index 80933bd44d2b9e9f1ff357d6e2969aa7e082abbb..ca51bf080045e0df06d706b60bc5fd1dca497d44 100644
--- a/tests/test_PyramidTransformation.cpp
+++ b/tests/test_PyramidTransformation.cpp
@@ -25,8 +25,7 @@ TEST_CASE("PyramidTransformation", "[geometry]")
const R3 d = {0, 3, 1};
const R3 e = {1, 2, 5};
- const std::array points = {a, b, c, d, e};
- const PyramidTransformation t(points);
+ const PyramidTransformation t(a, b, c, d, e);
SECTION("points")
{
diff --git a/tests/test_Q1Transformation.cpp b/tests/test_Q1Transformation.cpp
deleted file mode 100644
index caa46f5bf003ad503b620fb02fe6662cc027a27d..0000000000000000000000000000000000000000
--- a/tests/test_Q1Transformation.cpp
+++ /dev/null
@@ -1,238 +0,0 @@
-#include <catch2/catch_approx.hpp>
-#include <catch2/catch_test_macros.hpp>
-
-#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
-#include <analysis/QuadratureManager.hpp>
-#include <geometry/Q1Transformation.hpp>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("Q1Transformation", "[geometry]")
-{
- SECTION("1D")
- {
- using R1 = TinyVector<1>;
-
- const R1 a = 0.3;
- const R1 b = 2.7;
-
- const std::array points = {a, b};
-
- const Q1Transformation<1> t(points);
-
- 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(R1{-1}) == Catch::Approx(1.2));
- REQUIRE(t.jacobianDeterminant(R1{0}) == Catch::Approx(1.2));
- REQUIRE(t.jacobianDeterminant(R1{3}) == Catch::Approx(1.2));
- }
-
- SECTION("2D")
- {
- using R2 = TinyVector<2>;
-
- const R2 a_hat = {-1, -1};
- const R2 b_hat = {+1, -1};
- const R2 c_hat = {+1, +1};
- const R2 d_hat = {-1, +1};
-
- const R2 m_hat = zero;
-
- const R2 a = {0, 0};
- const R2 b = {8, -2};
- const R2 c = {12, 7};
- const R2 d = {3, 7};
-
- const R2 m = 0.25 * (a + b + c + d);
-
- const std::array points = {a, b, c, d};
- const Q1Transformation<2> t(points);
-
- SECTION("values")
- {
- REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
- REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
-
- REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
- REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
-
- REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
- REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
-
- REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
- REQUIRE(t(d_hat)[1] == Catch::Approx(d[1]));
-
- REQUIRE(t(m_hat)[0] == Catch::Approx(m[0]));
- REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
- }
-
- SECTION("Jacobian determinant")
- {
- SECTION("at points")
- {
- auto detJ = [](const R2& X) {
- const double x = X[0];
- const double y = X[1];
-
- return 0.25 * ((0.5 * x + 4) * (y + 17) - 0.5 * (x + 7) * (y - 1));
- };
-
- REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
- REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
- REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
- REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
-
- REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
- }
-
- SECTION("Gauss order 1")
- {
- // One point is enough in 2d
- const QuadratureFormula<2>& gauss =
- QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(1));
-
- double surface = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
- }
-
- // 71.5 is actually the surface of the quadrangle
- REQUIRE(surface == Catch::Approx(71.5));
- }
-
- SECTION("Gauss order 3")
- {
- auto p = [](const R2& X) {
- const double& x = X[0];
- const double& y = X[1];
-
- return (2 * x + 3 * y + 2) * (x - 2 * y - 1);
- };
-
- // Jacbian determinant is a degree 1 polynomial, so the
- // following formula is required to reach exactness
- const QuadratureFormula<2>& gauss =
- QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(3));
-
- double integral = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
- }
-
- REQUIRE(integral == Catch::Approx(-479. / 2));
- }
- }
- }
-
- SECTION("3D")
- {
- using R3 = TinyVector<3>;
-
- const R3 a_hat = {-1, -1, -1};
- const R3 b_hat = {+1, -1, -1};
- const R3 c_hat = {+1, +1, -1};
- const R3 d_hat = {-1, +1, -1};
- const R3 e_hat = {-1, -1, +1};
- const R3 f_hat = {+1, -1, +1};
- const R3 g_hat = {+1, +1, +1};
- const R3 h_hat = {-1, +1, +1};
-
- const R3 m_hat = zero;
-
- const R3 a = {0, 0, 0};
- const R3 b = {3, 1, 3};
- const R3 c = {2, 5, 2};
- const R3 d = {0, 3, 1};
- const R3 e = {1, 2, 5};
- const R3 f = {3, 1, 7};
- const R3 g = {2, 5, 5};
- const R3 h = {0, 3, 6};
-
- const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
-
- const std::array points = {a, b, c, d, e, f, g, h};
- const Q1Transformation<3> t(points);
-
- SECTION("values")
- {
- REQUIRE(l2Norm(t(a_hat) - a) == Catch::Approx(0));
- REQUIRE(l2Norm(t(b_hat) - b) == Catch::Approx(0));
- REQUIRE(l2Norm(t(c_hat) - c) == Catch::Approx(0));
- REQUIRE(l2Norm(t(d_hat) - d) == Catch::Approx(0));
- REQUIRE(l2Norm(t(e_hat) - e) == Catch::Approx(0));
- REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
- REQUIRE(l2Norm(t(g_hat) - g) == Catch::Approx(0));
- REQUIRE(l2Norm(t(h_hat) - h) == Catch::Approx(0));
-
- REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0));
- }
-
- SECTION("Jacobian determinant")
- {
- SECTION("at points")
- {
- auto detJ = [](const R3& X) {
- const double x = X[0];
- const double y = X[1];
- const double z = X[2];
-
- return -(3 * x * y * y * z + 9 * y * y * z - 6 * x * x * y * z + 6 * x * y * z - 96 * y * z + 14 * x * x * z -
- 49 * x * z + 119 * z - 21 * x * y * y + 21 * y * y + 90 * x * y - 10 * y + 40 * x * x - 109 * x -
- 491) /
- 128;
- };
-
- REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
- REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
- REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
- REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
- REQUIRE(t.jacobianDeterminant(e_hat) == Catch::Approx(detJ(e_hat)));
- REQUIRE(t.jacobianDeterminant(f_hat) == Catch::Approx(detJ(f_hat)));
- REQUIRE(t.jacobianDeterminant(g_hat) == Catch::Approx(detJ(g_hat)));
- REQUIRE(t.jacobianDeterminant(h_hat) == Catch::Approx(detJ(h_hat)));
-
- REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
- }
-
- SECTION("Gauss order 3")
- {
- // The jacobian determinant is of maximal degree 2 in each variable
- const QuadratureFormula<3>& gauss =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(2));
-
- double volume = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
- }
-
- // 353. / 12 is actually the volume of the hexahedron
- REQUIRE(volume == Catch::Approx(353. / 12));
- }
-
- SECTION("Gauss order 4")
- {
- 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 * y * z + z * z + 2 * x - 3 * y + 0.5 * z + 3;
- };
-
- // Jacbian determinant is a degree 2 polynomial, so the
- // following formula is required to reach exactness
- const QuadratureFormula<3>& gauss =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
-
- double integral = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
- }
-
- REQUIRE(integral == Catch::Approx(8795. / 72));
- }
- }
- }
-}
diff --git a/tests/test_SquareTransformation.cpp b/tests/test_SquareTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..19523959948a9ae3980e5df9f6efb4dcb0e12734
--- /dev/null
+++ b/tests/test_SquareTransformation.cpp
@@ -0,0 +1,104 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/SquareTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SquareTransformation", "[geometry]")
+{
+ using R2 = TinyVector<2>;
+
+ const R2 a_hat = {-1, -1};
+ const R2 b_hat = {+1, -1};
+ const R2 c_hat = {+1, +1};
+ const R2 d_hat = {-1, +1};
+
+ const R2 m_hat = zero;
+
+ const R2 a = {0, 0};
+ const R2 b = {8, -2};
+ const R2 c = {12, 7};
+ const R2 d = {3, 7};
+
+ const R2 m = 0.25 * (a + b + c + d);
+
+ const SquareTransformation t(a, b, c, d);
+
+ SECTION("values")
+ {
+ REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
+ REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
+
+ REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
+ REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
+
+ REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
+ REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
+
+ REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
+ REQUIRE(t(d_hat)[1] == Catch::Approx(d[1]));
+
+ REQUIRE(t(m_hat)[0] == Catch::Approx(m[0]));
+ REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
+ }
+
+ SECTION("Jacobian determinant")
+ {
+ SECTION("at points")
+ {
+ auto detJ = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+
+ return 0.25 * ((0.5 * x + 4) * (y + 17) - 0.5 * (x + 7) * (y - 1));
+ };
+
+ REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
+ REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
+ REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
+ REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
+
+ REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+ }
+
+ SECTION("Gauss order 1")
+ {
+ // One point is enough in 2d
+ const QuadratureFormula<2>& gauss =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(1));
+
+ double surface = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ }
+
+ // 71.5 is actually the surface of the quadrangle
+ REQUIRE(surface == Catch::Approx(71.5));
+ }
+
+ SECTION("Gauss order 3")
+ {
+ auto p = [](const R2& X) {
+ const double& x = X[0];
+ const double& y = X[1];
+
+ return (2 * x + 3 * y + 2) * (x - 2 * y - 1);
+ };
+
+ // Jacbian determinant is a degree 1 polynomial, so the
+ // following formula is required to reach exactness
+ const QuadratureFormula<2>& gauss =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(3));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ }
+
+ REQUIRE(integral == Catch::Approx(-479. / 2));
+ }
+ }
+}
diff --git a/tests/test_TetrahedronGaussQuadrature.cpp b/tests/test_TetrahedronGaussQuadrature.cpp
index 00f3eecf101e6e9bbc45ba3866cc06237f922569..605d1dc1fa7476719fe16ac595ff55491f5bd8e8 100644
--- a/tests/test_TetrahedronGaussQuadrature.cpp
+++ b/tests/test_TetrahedronGaussQuadrature.cpp
@@ -5,7 +5,7 @@
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <analysis/TetrahedronGaussQuadrature.hpp>
-#include <geometry/AffineTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
@@ -20,7 +20,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
const R3 C{-1, +1, -1};
const R3 D{-1, -1, +1};
- AffineTransformation<3> t{{A, B, C, D}};
+ TetrahedronTransformation t{A, B, C, D};
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -33,8 +33,9 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
return t.jacobianDeterminant() * value;
};
- auto integrate_on_tetra = [](auto f, auto quadrature_formula, const std::array<TinyVector<3>, 4>& tetrahedron) {
- AffineTransformation<3> t{tetrahedron};
+ auto integrate_on_tetra = [](auto f, auto quadrature_formula, const TinyVector<3>& a, const TinyVector<3>& b,
+ const TinyVector<3>& c, const TinyVector<3>& d) {
+ TetrahedronTransformation t{a, b, c, d};
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -64,14 +65,14 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
const double int_T_hat = integrate(f, quadrature_formula);
const double int_refined //
- = integrate_on_tetra(f, quadrature_formula, {A, M_AB, M_AC, M_AD}) +
- integrate_on_tetra(f, quadrature_formula, {B, M_AB, M_BD, M_BC}) +
- integrate_on_tetra(f, quadrature_formula, {C, M_AC, M_BC, M_CD}) +
- integrate_on_tetra(f, quadrature_formula, {D, M_AD, M_CD, M_BD}) +
- integrate_on_tetra(f, quadrature_formula, {M_BD, M_AC, M_AD, M_CD}) +
- integrate_on_tetra(f, quadrature_formula, {M_AB, M_AC, M_AD, M_BD}) +
- integrate_on_tetra(f, quadrature_formula, {M_BC, M_AB, M_BD, M_AC}) +
- integrate_on_tetra(f, quadrature_formula, {M_BC, M_AC, M_BD, M_CD});
+ = integrate_on_tetra(f, quadrature_formula, A, M_AB, M_AC, M_AD) +
+ integrate_on_tetra(f, quadrature_formula, B, M_AB, M_BD, M_BC) +
+ integrate_on_tetra(f, quadrature_formula, C, M_AC, M_BC, M_CD) +
+ integrate_on_tetra(f, quadrature_formula, D, M_AD, M_CD, M_BD) +
+ integrate_on_tetra(f, quadrature_formula, M_BD, M_AC, M_AD, M_CD) +
+ integrate_on_tetra(f, quadrature_formula, M_AB, M_AC, M_AD, M_BD) +
+ integrate_on_tetra(f, quadrature_formula, M_BC, M_AB, M_BD, M_AC) +
+ integrate_on_tetra(f, quadrature_formula, M_BC, M_AC, M_BD, M_CD);
return -std::log(std::abs(int_refined - exact_value) / std::abs(int_T_hat - exact_value)) / std::log(2);
};
diff --git a/tests/test_TetrahedronTransformation.cpp b/tests/test_TetrahedronTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a9d20275fd1df42db4e5ea723b585c2872e69c4
--- /dev/null
+++ b/tests/test_TetrahedronTransformation.cpp
@@ -0,0 +1,40 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <geometry/TetrahedronTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("TetrahedronTransformation", "[geometry]")
+{
+ using R3 = TinyVector<3>;
+
+ const R3 a = {1, 2, 1};
+ const R3 b = {3, 1, 3};
+ const R3 c = {2, 5, 2};
+ const R3 d = {2, 3, 4};
+
+ const TetrahedronTransformation t(a, b, c, d);
+
+ REQUIRE(t({0, 0, 0})[0] == Catch::Approx(1));
+ REQUIRE(t({0, 0, 0})[1] == Catch::Approx(2));
+ REQUIRE(t({0, 0, 0})[2] == Catch::Approx(1));
+
+ REQUIRE(t({1, 0, 0})[0] == Catch::Approx(3));
+ REQUIRE(t({1, 0, 0})[1] == Catch::Approx(1));
+ REQUIRE(t({1, 0, 0})[2] == Catch::Approx(3));
+
+ REQUIRE(t({0, 1, 0})[0] == Catch::Approx(2));
+ REQUIRE(t({0, 1, 0})[1] == Catch::Approx(5));
+ REQUIRE(t({0, 1, 0})[2] == Catch::Approx(2));
+
+ REQUIRE(t({0, 0, 1})[0] == Catch::Approx(2));
+ REQUIRE(t({0, 0, 1})[1] == Catch::Approx(3));
+ REQUIRE(t({0, 0, 1})[2] == Catch::Approx(4));
+
+ REQUIRE(t({0.25, 0.25, 0.25})[0] == Catch::Approx(2));
+ REQUIRE(t({0.25, 0.25, 0.25})[1] == Catch::Approx(11. / 4));
+ REQUIRE(t({0.25, 0.25, 0.25})[2] == Catch::Approx(2.5));
+
+ REQUIRE(t.jacobianDeterminant() == Catch::Approx(14));
+}
diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp
index 879c47b94b3326e96bbefecc40047ea50fb49b71..0bb7f34ec3cdc057249f9fa2c0f3d9eaf6dc81c5 100644
--- a/tests/test_TriangleGaussQuadrature.cpp
+++ b/tests/test_TriangleGaussQuadrature.cpp
@@ -5,7 +5,7 @@
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <analysis/TriangleGaussQuadrature.hpp>
-#include <geometry/AffineTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
@@ -19,7 +19,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
const R2 B{+1, -1};
const R2 C{-1, +1};
- AffineTransformation<2> t{{A, B, C}};
+ TriangleTransformation t{A, B, C};
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -32,8 +32,9 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
return t.jacobianDeterminant() * value;
};
- auto integrate_on_triangle = [](auto f, auto quadrature_formula, const std::array<TinyVector<2>, 3>& triangle) {
- AffineTransformation<2> t{triangle};
+ 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);
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -50,10 +51,10 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
using R2 = TinyVector<2>;
const double int_T_hat = integrate(f, quadrature_formula);
const double int_refined //
- = integrate_on_triangle(f, quadrature_formula, {R2{-1, -1}, R2{0, -1}, R2{-1, 0}}) +
- integrate_on_triangle(f, quadrature_formula, {R2{0, -1}, R2{0, 0}, R2{-1, 0}}) +
- integrate_on_triangle(f, quadrature_formula, {R2{0, -1}, R2{1, -1}, R2{0, 0}}) +
- integrate_on_triangle(f, quadrature_formula, {R2{-1, 0}, R2{0, 0}, R2{-1, 1}});
+ = integrate_on_triangle(f, quadrature_formula, R2{-1, -1}, R2{0, -1}, R2{-1, 0}) +
+ integrate_on_triangle(f, quadrature_formula, R2{0, -1}, R2{0, 0}, R2{-1, 0}) +
+ integrate_on_triangle(f, quadrature_formula, R2{0, -1}, R2{1, -1}, R2{0, 0}) +
+ integrate_on_triangle(f, quadrature_formula, R2{-1, 0}, R2{0, 0}, R2{-1, 1});
return -std::log(std::abs(int_refined - exact_value) / std::abs(int_T_hat - exact_value)) / std::log(2);
};
diff --git a/tests/test_TriangleTransformation.cpp b/tests/test_TriangleTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1c05aabc7be130e9b36bc4b52a3ebe0733e8456
--- /dev/null
+++ b/tests/test_TriangleTransformation.cpp
@@ -0,0 +1,31 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <geometry/TriangleTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("TriangleTransformation", "[geometry]")
+{
+ using R2 = TinyVector<2>;
+
+ const R2 a = {1, 2};
+ const R2 b = {3, 1};
+ const R2 c = {2, 5};
+
+ const TriangleTransformation t(a, b, c);
+
+ 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({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.jacobianDeterminant() == Catch::Approx(7));
+}