diff --git a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
index caf0762e9fb297af55eef90269a4b5d659a97b9a..4c3748fd173650dade28e5759b9b078cc73ff126 100644
--- a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
+++ b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
@@ -178,9 +178,8 @@ void rotateGeometryToXY(std::vector<GeoLib::Point*>& points,
auto const [plane_normal, d] =
GeoLib::getNewellPlane(points.begin(), points.end());
// rotate points into x-y-plane
- // Todo (TF) Remove when rotateGeometryToXY uses Eigen for rot_mat
- Eigen::Matrix3d rotation_matrix_;
- GeoLib::computeRotationMatrixToXY(plane_normal, rotation_matrix_);
+ Eigen::Matrix3d const rotation_matrix_ =
+ GeoLib::computeRotationMatrixToXY(plane_normal);
GeoLib::rotatePoints(rotation_matrix_, points.begin(), points.end());
// Todo (TF) Remove when rotateGeometryToXY uses Eigen for rot_mat
for (int r = 0; r < 3; r++)
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index 39a0b96666ff32b24b9dac66e8a019470aadbd81..76090a3675deb1ee9ad0b98cc659fe18215e14c0 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -117,80 +117,6 @@ void compute3DRotationMatrixToX(MathLib::Vector3 const& v,
}
}
-template <class T_MATRIX>
-void computeRotationMatrixToXY(Eigen::Vector3d const& n, T_MATRIX& rot_mat)
-{
- // check if normal points already in the right direction
- if (n[0] == 0 && n[1] == 0) {
- rot_mat(0,1) = 0.0;
- rot_mat(0,2) = 0.0;
- rot_mat(1,0) = 0.0;
- rot_mat(1,1) = 1.0;
- rot_mat(1,2) = 0.0;
- rot_mat(2,0) = 0.0;
- rot_mat(2,1) = 0.0;
-
- if (n[2] > 0) {
- // identity matrix
- rot_mat(0,0) = 1.0;
- rot_mat(2,2) = 1.0;
- } else {
- // rotate by pi about the y-axis
- rot_mat(0,0) = -1.0;
- rot_mat(2,2) = -1.0;
- }
-
- return;
- }
-
- // sqrt (n_1^2 + n_3^2)
- double const h0(sqrt(n[0]*n[0]+n[2]*n[2]));
-
- // In case the x and z components of the normal are both zero the rotation
- // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
- // required. The angle is either pi/2 or 3/2*pi. Thus the components of
- // rot_mat are as follows.
- if (h0 < std::numeric_limits<double>::epsilon()) {
- rot_mat(0,0) = 1.0;
- rot_mat(0,1) = 0.0;
- rot_mat(0,2) = 0.0;
- rot_mat(1,0) = 0.0;
- rot_mat(1,1) = 0.0;
- if (n[1] > 0)
- {
- rot_mat(1,2) = -1.0;
- }
- else
- {
- rot_mat(1, 2) = 1.0;
- }
- rot_mat(2,0) = 0.0;
- if (n[1] > 0)
- {
- rot_mat(2,1) = 1.0;
- }
- else
- {
- rot_mat(2, 1) = -1.0;
- }
- rot_mat(2,2) = 0.0;
- return;
- }
-
- double const h1(1 / n.norm());
-
- // general case: calculate entries of rotation matrix
- rot_mat(0, 0) = n[2] / h0;
- rot_mat(0, 1) = 0;
- rot_mat(0, 2) = - n[0] / h0;
- rot_mat(1, 0) = - n[1]*n[0]/h0 * h1;
- rot_mat(1, 1) = h0 * h1;
- rot_mat(1, 2) = - n[1]*n[2]/h0 * h1;
- rot_mat(2, 0) = n[0] * h1;
- rot_mat(2, 1) = n[1] * h1;
- rot_mat(2, 2) = n[2] * h1;
-}
-
template <typename InputIterator>
void rotatePoints(Eigen::Matrix3d const& rot_mat, InputIterator pnts_begin,
InputIterator pnts_end)
@@ -215,8 +141,7 @@ Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin,
GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
// rotate points into x-y-plane
- Eigen::Matrix3d rot_mat;
- computeRotationMatrixToXY(plane_normal, rot_mat);
+ Eigen::Matrix3d const rot_mat = computeRotationMatrixToXY(plane_normal);
rotatePoints(rot_mat, r_pnts_begin, r_pnts_end);
for (auto it = r_pnts_begin; it != r_pnts_end; ++it)
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 4feae3e931a4c0ac6c573a80664369236f14f5b4..509ec76ba0374bee16f961775a62755216e0818f 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -292,6 +292,69 @@ void rotatePoints(Eigen::Matrix3d const& rot_mat,
rotatePoints(rot_mat, pnts.begin(), pnts.end());
}
+Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const& n)
+{
+ Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
+ // check if normal points already in the right direction
+ if (n[0] == 0 && n[1] == 0)
+ {
+ rot_mat(1, 1) = 1.0;
+
+ if (n[2] > 0)
+ {
+ // identity matrix
+ rot_mat(0, 0) = 1.0;
+ rot_mat(2, 2) = 1.0;
+ }
+ else
+ {
+ // rotate by pi about the y-axis
+ rot_mat(0, 0) = -1.0;
+ rot_mat(2, 2) = -1.0;
+ }
+
+ return rot_mat;
+ }
+
+ // sqrt (n_1^2 + n_3^2)
+ double const h0(std::sqrt(n[0] * n[0] + n[2] * n[2]));
+
+ // In case the x and z components of the normal are both zero the rotation
+ // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
+ // required. The angle is either pi/2 or 3/2*pi. Thus the components of
+ // rot_mat are as follows.
+ if (h0 < std::numeric_limits<double>::epsilon())
+ {
+ rot_mat(0, 0) = 1.0;
+ if (n[1] > 0)
+ {
+ rot_mat(1, 2) = -1.0;
+ rot_mat(2, 1) = 1.0;
+ }
+ else
+ {
+ rot_mat(1, 2) = 1.0;
+ rot_mat(2, 1) = -1.0;
+ }
+ return rot_mat;
+ }
+
+ double const h1(1 / n.norm());
+
+ // general case: calculate entries of rotation matrix
+ rot_mat(0, 0) = n[2] / h0;
+ rot_mat(0, 1) = 0;
+ rot_mat(0, 2) = -n[0] / h0;
+ rot_mat(1, 0) = -n[1] * n[0] / h0 * h1;
+ rot_mat(1, 1) = h0 * h1;
+ rot_mat(1, 2) = -n[1] * n[2] / h0 * h1;
+ rot_mat(2, 0) = n[0] * h1;
+ rot_mat(2, 1) = n[1] * h1;
+ rot_mat(2, 2) = n[2] * h1;
+
+ return rot_mat;
+}
+
Eigen::Matrix3d rotatePointsToXY(std::vector<GeoLib::Point*>& pnts)
{
return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
@@ -382,8 +445,8 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
// 2 rotate points
double d_polygon;
std::tie(plane_normal, d_polygon) = GeoLib::getNewellPlane(*polygon_pnts);
- Eigen::Matrix3d rot_mat;
- GeoLib::computeRotationMatrixToXY(plane_normal, rot_mat);
+ Eigen::Matrix3d const rot_mat =
+ GeoLib::computeRotationMatrixToXY(plane_normal);
GeoLib::rotatePoints(rot_mat, *polygon_pnts);
// 3 set z coord to zero
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index cb3f4b8c19c1c6a86b31caa8c94aa457d6074e2d..5de5f0fa5cd1cf1beeb17f1619043a169c4ed938 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -100,10 +100,9 @@ void compute3DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX& rot_mat);
* Method computes the rotation matrix that rotates the given vector parallel to
* the \f$z\f$ axis.
* @param n the (3d) vector that is rotated parallel to the \f$z\f$ axis
- * @param rot_mat 3x3 rotation matrix
+ * @return rot_mat 3x3 rotation matrix
*/
-template <class T_MATRIX>
-void computeRotationMatrixToXY(Eigen::Vector3d const& n, T_MATRIX& rot_mat);
+Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const& n);
/**
* rotate points according to the rotation matrix
diff --git a/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
index b8ece6590dc3114e8ba70382794fa0ef4021df14..33dd3cd879e3bb2bf787387223baaeae414b808a 100644
--- a/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
+++ b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
@@ -16,9 +16,6 @@
#include "GeoLib/AnalyticalGeometry.h"
#include "GeoLib/Polygon.h"
-#include "MathLib/Vector3.h"
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-
#include "MeshLib/Node.h"
namespace MeshGeoToolsLib
@@ -38,8 +35,7 @@ std::vector<bool> markNodesOutSideOfPolygon(
rotated_nodes.push_back(new GeoLib::Point(*node, node->getID()));
}
// 2 rotate the Points
- Eigen::Matrix3d rot_mat;
- GeoLib::computeRotationMatrixToXY(normal, rot_mat);
+ Eigen::Matrix3d const rot_mat = GeoLib::computeRotationMatrixToXY(normal);
GeoLib::rotatePoints(rot_mat, rotated_nodes);
// 3 set z coord to zero
std::for_each(rotated_nodes.begin(), rotated_nodes.end(),
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index 9f683adcaf1159a04278a34bbe35533aa14bfc96..1d548a41d207fb5a4b1bd304c0e99b3cf46e630a 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -59,7 +59,7 @@ void getRotationMatrixToGlobal(const unsigned element_dimension,
// get plane normal
auto const [plane_normal, d] = GeoLib::getNewellPlane(points);
// compute a rotation matrix to XY
- GeoLib::computeRotationMatrixToXY(plane_normal, matR);
+ matR = GeoLib::computeRotationMatrixToXY(plane_normal);
// set a transposed matrix
matR.transposeInPlace();
}
diff --git a/Tests/GeoLib/TestComputeRotationMatrix.cpp b/Tests/GeoLib/TestComputeRotationMatrix.cpp
index b2fd851ecf18e7cfaad241d250c3dbe05ecd650b..092f886c386bfc87e71fbe2603f7074f8e401816 100644
--- a/Tests/GeoLib/TestComputeRotationMatrix.cpp
+++ b/Tests/GeoLib/TestComputeRotationMatrix.cpp
@@ -13,21 +13,18 @@
#include "GeoLib/AnalyticalGeometry.h"
-auto test3equal = [](double a, double b, double c, double const* result)
-{
+auto test3equal = [](double a, double b, double c,
+ Eigen::Vector3d const& result) {
EXPECT_EQ(a, result[0]);
EXPECT_EQ(b, result[1]);
EXPECT_EQ(c, result[2]);
-
- delete[] result;
};
TEST(GeoLib, ComputeRotationMatrixToXYnegative)
{
Eigen::Vector3d const n({0.0, -1.0, 0.0});
- MathLib::DenseMatrix<double> rot_mat(3,3,0.0);
+ Eigen::Matrix3d const rot_mat = GeoLib::computeRotationMatrixToXY(n);
- GeoLib::computeRotationMatrixToXY(n, rot_mat);
EXPECT_EQ(1.0, rot_mat(0,0));
EXPECT_EQ(0.0, rot_mat(0,1));
EXPECT_EQ(0.0, rot_mat(0,2));
@@ -38,22 +35,21 @@ TEST(GeoLib, ComputeRotationMatrixToXYnegative)
EXPECT_EQ(-1.0, rot_mat(2,1));
EXPECT_EQ(0.0, rot_mat(2,2));
- MathLib::Vector3 const x(0.0,1.0,0.0);
- test3equal(0, 0, -1, rot_mat*x.getCoords());
+ Eigen::Vector3d const x({0.0, 1.0, 0.0});
+ test3equal(0, 0, -1, rot_mat * x);
- MathLib::Vector3 const x0(10.0,1.0,0.0);
- test3equal(10, 0, -1, rot_mat*x0.getCoords());
+ Eigen::Vector3d const x0({10.0, 1.0, 0.0});
+ test3equal(10, 0, -1, rot_mat * x0);
- MathLib::Vector3 const x1(10.0,0.0,10.0);
- test3equal(10, 10, 0, rot_mat*x1.getCoords());
+ Eigen::Vector3d const x1({10.0, 0.0, 10.0});
+ test3equal(10, 10, 0, rot_mat * x1);
}
TEST(GeoLib, ComputeRotationMatrixToXYpositive)
{
Eigen::Vector3d const n{0.0, 1.0, 0.0};
- MathLib::DenseMatrix<double> rot_mat(3,3,0.0);
+ Eigen::Matrix3d const rot_mat = GeoLib::computeRotationMatrixToXY(n);
- GeoLib::computeRotationMatrixToXY(n, rot_mat);
EXPECT_EQ(1.0, rot_mat(0,0));
EXPECT_EQ(0.0, rot_mat(0,1));
EXPECT_EQ(0.0, rot_mat(0,2));
@@ -64,13 +60,13 @@ TEST(GeoLib, ComputeRotationMatrixToXYpositive)
EXPECT_EQ(1.0, rot_mat(2,1));
EXPECT_EQ(0.0, rot_mat(2,2));
- MathLib::Vector3 const x(0.0,1.0,0.0);
- test3equal(0, 0, 1, rot_mat*x.getCoords());
+ Eigen::Vector3d const x(0.0, 1.0, 0.0);
+ test3equal(0, 0, 1, rot_mat * x);
- MathLib::Vector3 const x0(10.0,1.0,0.0);
- test3equal(10, 0, 1, rot_mat*x0.getCoords());
+ Eigen::Vector3d const x0(10.0, 1.0, 0.0);
+ test3equal(10, 0, 1, rot_mat * x0);
- MathLib::Vector3 const x1(10.0,0.0,10.0);
- test3equal(10, -10, 0, rot_mat*x1.getCoords());
+ Eigen::Vector3d const x1(10.0, 0.0, 10.0);
+ test3equal(10, -10, 0, rot_mat * x1);
}