From 891fea9676e4e64b6fad362f78687f248bfc0b6c Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 10 Dec 2020 10:59:39 +0100
Subject: [PATCH] [GL] Use Eigen::Vector3d in getNewellPlane impl.
---
.../MeshGeoTools/VerticalSliceFromLayers.cpp | 5 +++-
GeoLib/AnalyticalGeometry-impl.h | 29 ++++++++++---------
GeoLib/AnalyticalGeometry.cpp | 8 +++--
GeoLib/AnalyticalGeometry.h | 10 +++----
MeshLib/ElementCoordinatesMappingLocal.cpp | 9 +++---
5 files changed, 36 insertions(+), 25 deletions(-)
diff --git a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
index e6a2d691783..e7cf307cfe7 100644
--- a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
+++ b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
@@ -196,8 +196,11 @@ void rotateGeometryToXY(std::vector<GeoLib::Point*>& points,
// compute the plane normal
auto const [plane_normal, d] =
GeoLib::getNewellPlane(points.begin(), points.end());
+ // ToDo (TF) remove when computeRotationMatrixToXY uses Eigen
+ MathLib::Vector3 const tmp_plane_normal(
+ {plane_normal[0], plane_normal[1], plane_normal[2]});
// rotate points into x-y-plane
- GeoLib::computeRotationMatrixToXY(plane_normal, rotation_matrix);
+ GeoLib::computeRotationMatrixToXY(tmp_plane_normal, rotation_matrix);
GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
GeoLib::AABB aabb(points.begin(), points.end());
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index 12086e2cf98..9944d7e0035 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -11,11 +11,11 @@
namespace GeoLib
{
template <typename InputIterator>
-std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
- InputIterator pnts_end)
+std::pair<Eigen::Vector3d, double> getNewellPlane(InputIterator pnts_begin,
+ InputIterator pnts_end)
{
- MathLib::Vector3 plane_normal;
- MathLib::Vector3 centroid;
+ Eigen::Vector3d plane_normal({0, 0, 0});
+ Eigen::Vector3d centroid({0, 0, 0});
for (auto i=std::prev(pnts_end), j=pnts_begin; j!=pnts_end; i = j, ++j) {
auto &pt_i = *(*i);
auto &pt_j = *(*j);
@@ -26,7 +26,7 @@ std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
plane_normal[2] += (pt_i[0] - pt_j[0])
* (pt_i[1] + pt_j[1]); // projection on xy
- centroid += MathLib::Vector3(pt_j);
+ centroid += Eigen::Map<Eigen::Vector3d const>(pt_j.getCoords());
}
plane_normal.normalize();
@@ -35,24 +35,24 @@ std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
double d = 0.0;
if (n_pnts > 0)
{
- d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts;
+ d = centroid.dot(plane_normal) / n_pnts;
}
return std::make_pair(plane_normal, d);
}
template <class T_POINT>
-std::pair<MathLib::Vector3, double> getNewellPlane(
+std::pair<Eigen::Vector3d, double> getNewellPlane(
const std::vector<T_POINT*>& pnts)
{
return getNewellPlane(pnts.begin(), pnts.end());
}
template <class T_POINT>
-std::pair<MathLib::Vector3, double> getNewellPlane(
+std::pair<Eigen::Vector3d, double> getNewellPlane(
const std::vector<T_POINT>& pnts)
{
- MathLib::Vector3 plane_normal;
- MathLib::Vector3 centroid;
+ Eigen::Vector3d plane_normal({0, 0, 0});
+ Eigen::Vector3d centroid({0, 0, 0});
std::size_t n_pnts(pnts.size());
for (std::size_t i = n_pnts - 1, j = 0; j < n_pnts; i = j, j++) {
plane_normal[0] += (pnts[i][1] - pnts[j][1])
@@ -62,11 +62,11 @@ std::pair<MathLib::Vector3, double> getNewellPlane(
plane_normal[2] += (pnts[i][0] - pnts[j][0])
* (pnts[i][1] + pnts[j][1]); // projection on xy
- centroid += MathLib::Vector3(pnts[j]);
+ centroid += Eigen::Map<Eigen::Vector3d const>(pnts[j].getCoords());
}
plane_normal.normalize();
- double d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts;
+ double const d = centroid.dot(plane_normal) / n_pnts;
return std::make_pair(plane_normal, d);
}
@@ -219,10 +219,13 @@ MathLib::DenseMatrix<double> rotatePointsToXY(InputIterator1 p_pnts_begin,
// compute the plane normal
auto const [plane_normal, d] =
GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
+ // ToDo (TF) remove when computeRotationMatrixToXY uses Eigen
+ MathLib::Vector3 const tmp_plane_normal(
+ {plane_normal[0], plane_normal[1], plane_normal[2]});
// rotate points into x-y-plane
MathLib::DenseMatrix<double> rot_mat(3, 3);
- computeRotationMatrixToXY(plane_normal, rot_mat);
+ computeRotationMatrixToXY(tmp_plane_normal, rot_mat);
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 9ee2d4b28bd..198c0443744 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -373,8 +373,12 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
}
// 2 rotate points
- double d_polygon (0.0);
- std::tie(plane_normal, d_polygon) = GeoLib::getNewellPlane(*polygon_pnts);
+ auto const [normal, d_polygon] = GeoLib::getNewellPlane(*polygon_pnts);
+ // ToDo (TF) remove when computeRotationMatrixToXY uses Eigen
+ plane_normal[0] = normal[0];
+ plane_normal[1] = normal[1];
+ plane_normal[2] = normal[2];
+
MathLib::DenseMatrix<double> rot_mat(3,3);
GeoLib::computeRotationMatrixToXY(plane_normal, rot_mat);
GeoLib::rotatePoints(rot_mat, *polygon_pnts);
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index d091d1e5307..b1a9d2ea727 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -58,8 +58,8 @@ Orientation getOrientationFast(MathLib::Point3d const& p0,
* polygon is located in, d parameter from the plane equation
*/
template <typename InputIterator>
-std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
- InputIterator pnts_end);
+std::pair<Eigen::Vector3d, double> getNewellPlane(InputIterator pnts_begin,
+ InputIterator pnts_end);
/**
* compute a supporting plane (represented by plane_normal and the value d) for the polygon
@@ -71,13 +71,13 @@ std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
* polygon is located in, parameter d from the plane equation
*/
template <class T_POINT>
-std::pair<MathLib::Vector3, double> getNewellPlane(
+std::pair<Eigen::Vector3d, double> getNewellPlane(
const std::vector<T_POINT*>& pnts);
-/** Same as getNewellPlane(pnts, plane_normal, d).
+/** Same as getNewellPlane(pnts).
*/
template <class T_POINT>
-std::pair<MathLib::Vector3, double> getNewellPlane(
+std::pair<Eigen::Vector3d, double> getNewellPlane(
const std::vector<T_POINT>& pnts);
/**
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index f4ac3b9ba7b..9b1b9151e24 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -57,12 +57,13 @@ void getRotationMatrixToGlobal(const unsigned element_dimension,
else if (global_dim == 3 && element_dimension == 2)
{
// get plane normal
- MathLib::Vector3 plane_normal;
- double d;
- std::tie(plane_normal, d) = GeoLib::getNewellPlane(points);
+ auto const [plane_normal, d] = GeoLib::getNewellPlane(points);
+ // ToDo (TF) remove when computeRotationMatrixToXY uses Eigen
+ MathLib::Vector3 const tmp_plane_normal(
+ {plane_normal[0], plane_normal[1], plane_normal[2]});
// compute a rotation matrix to XY
- GeoLib::computeRotationMatrixToXY(plane_normal, matR);
+ GeoLib::computeRotationMatrixToXY(tmp_plane_normal, matR);
// set a transposed matrix
matR.transposeInPlace();
}
--
GitLab