From d47161e759a29804ab8d9bcd25c069dd7b9f0115 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 21 Feb 2020 17:02:34 +0100
Subject: [PATCH] [ParL] Add transformation_3d to coordinate system.

This is explicitly allowing the 2d coordinate system
being given (contrary to the transformation<3> specialization).
---
 ParameterLib/CoordinateSystem.cpp | 25 +++++++++++++++++++++++++
 ParameterLib/CoordinateSystem.h   |  3 +++
 2 files changed, 28 insertions(+)

diff --git a/ParameterLib/CoordinateSystem.cpp b/ParameterLib/CoordinateSystem.cpp
index 837e6659604..0dc9fc8bb4a 100644
--- a/ParameterLib/CoordinateSystem.cpp
+++ b/ParameterLib/CoordinateSystem.cpp
@@ -126,6 +126,31 @@ Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation<3>(
     return t;
 }
 
+Eigen::Matrix<double, 3, 3> CoordinateSystem::transformation_3d(
+    SpatialPosition const& pos) const
+{
+    if (_base[2] != nullptr)
+    {
+        return transformation<3>(pos);
+    }
+
+    auto e0 = (*_base[0])(0 /* time independent */, pos);
+    auto e1 = (*_base[1])(0 /* time independent */, pos);
+    Eigen::Matrix<double, 3, 3> t = Eigen::Matrix<double, 3, 3>::Identity();
+    t.template topLeftCorner<2, 2>() << e0[0], e1[0], e0[1], e1[1];
+
+#ifndef NDEBUG
+    if (std::abs(t.determinant() - 1) > std::numeric_limits<double>::epsilon())
+    {
+        OGS_FATAL(
+            "The determinant of the coordinate system transformation matrix is "
+            "'%g', which is not sufficiently close to unity.",
+            t.determinant());
+    }
+#endif  // NDEBUG
+    return t;
+}
+
 template <int Dimension>
 Eigen::Matrix<double, Dimension, Dimension> CoordinateSystem::rotateTensor(
     std::vector<double> const& values, SpatialPosition const& pos) const
diff --git a/ParameterLib/CoordinateSystem.h b/ParameterLib/CoordinateSystem.h
index 860647375e4..bc679004168 100644
--- a/ParameterLib/CoordinateSystem.h
+++ b/ParameterLib/CoordinateSystem.h
@@ -35,6 +35,9 @@ struct CoordinateSystem final
     Eigen::Matrix<double, Dimension, Dimension> transformation(
         SpatialPosition const& pos) const;
 
+    Eigen::Matrix<double, 3, 3> transformation_3d(
+        SpatialPosition const& pos) const;
+
     template <int Dimension>
     Eigen::Matrix<double, Dimension, Dimension> rotateTensor(
         std::vector<double> const& values, SpatialPosition const& pos) const;
-- 
GitLab