diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index d867f34ca5225cc2cc3f87c42d60d4411eabf2e7..a646b6f40dd9dfe5e4aea7b1ee4f3e609c74ed4e 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -70,59 +70,6 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(
return std::make_pair(plane_normal, d);
}
-template <class T_MATRIX>
-void compute2DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat)
-{
- const double cos_theta = v[0];
- const double sin_theta = v[1];
- rot_mat(0,0) = rot_mat(1,1) = cos_theta;
- rot_mat(0,1) = sin_theta;
- rot_mat(1,0) = -sin_theta;
- rot_mat(0,2) = rot_mat(1,2) = rot_mat(2,0) = rot_mat(2,1) = 0.0;
- rot_mat(2,2) = 1.0;
-}
-
-template <class T_MATRIX>
-void compute3DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat)
-{
- double const eps = std::numeric_limits<double>::epsilon();
- // a vector on the plane
- Eigen::Vector3d yy({0, 0, 0});
- if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
- {
- yy[2] = 1.0;
- }
- else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
- {
- yy[0] = 1.0;
- }
- else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
- {
- yy[1] = 1.0;
- }
- else
- {
- for (unsigned i = 0; i < 3; i++)
- {
- if (std::abs(v[i]) > 0.0)
- {
- yy[i] = -v[i];
- break;
- }
- }
- }
- // z"_vec
- Eigen::Vector3d const zz = v.cross(yy).normalized();
- // y"_vec
- Eigen::Vector3d const y = zz.cross(v).normalized();
-
- for (unsigned i = 0; i < 3; ++i)
- {
- rot_mat(0, i) = v[i];
- rot_mat(1, i) = y[i];
- rot_mat(2, i) = zz[i];
- }
-}
template <typename InputIterator>
void rotatePoints(Eigen::Matrix3d const& rot_mat, InputIterator pnts_begin,
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index d59a9fc20767f4848222bce41c86dfbe54cec695..8d1b5d3530d860e5d713429ffa9565f9aa347927 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -681,4 +681,56 @@ void sortSegments(
}
}
+Eigen::Matrix3d compute2DRotationMatrixToX(Eigen::Vector3d const& v)
+{
+ Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
+ const double cos_theta = v[0];
+ const double sin_theta = v[1];
+ rot_mat(0,0) = rot_mat(1,1) = cos_theta;
+ rot_mat(0,1) = sin_theta;
+ rot_mat(1,0) = -sin_theta;
+ rot_mat(2,2) = 1.0;
+ return rot_mat;
+}
+
+Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const& v)
+{
+ // a vector on the plane
+ Eigen::Vector3d yy = Eigen::Vector3d::Zero();
+ auto const eps = std::numeric_limits<double>::epsilon();
+ if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
+ {
+ yy[2] = 1.0;
+ }
+ else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
+ {
+ yy[0] = 1.0;
+ }
+ else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
+ {
+ yy[1] = 1.0;
+ }
+ else
+ {
+ for (unsigned i = 0; i < 3; i++)
+ {
+ if (std::abs(v[i]) > 0.0)
+ {
+ yy[i] = -v[i];
+ break;
+ }
+ }
+ }
+ // z"_vec
+ Eigen::Vector3d const zz = v.cross(yy).normalized();
+ // y"_vec
+ yy = zz.cross(v).normalized();
+
+ Eigen::Matrix3d rot_mat;
+ rot_mat.row(0) = v;
+ rot_mat.row(1) = yy;
+ rot_mat.row(2) = zz;
+ return rot_mat;
+}
+
} // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index f911cedf91f6a0c1e3c9b340e8cc25f477e52c07..08f8ca68d74a9676613c0e9769a18242fa812b84 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -78,20 +78,20 @@ std::pair<Eigen::Vector3d, double> getNewellPlane(
const std::vector<T_POINT>& pnts);
/**
- * Computes a rotation matrix that rotates the given 2D normal vector parallel to X-axis
- * @param v a 2D normal vector to be rotated
- * @param rot_mat a 2x2 rotation matrix
+ * Computes a rotation matrix that rotates the given 2D normal vector parallel
+ * to X-axis
+ * @param v a 2D normal vector to be rotated
+ * @return a 3x3 rotation matrix where the upper, left, 2x2 block
+ * contains the entries necessary for the 2D rotation
*/
-template<class T_MATRIX>
-void compute2DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX & rot_mat);
+Eigen::Matrix3d compute2DRotationMatrixToX(Eigen::Vector3d const& v);
/**
* Computes a rotation matrix that rotates the given 3D normal vector parallel to X-axis
* @param v a 3D normal vector to be rotated
- * @param rot_mat a 3x3 rotation matrix
+ * @return a 3x3 rotation matrix
*/
-template <class T_MATRIX>
-void compute3DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat);
+Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const& v);
/**
* Method computes the rotation matrix that rotates the given vector parallel to
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index bcd654bb25d92cd7a2f5ffe7fd974a1a87bdcd13..730cdfe732dddf4c3763a0ef857e3fc5e4fc4459 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -32,26 +32,27 @@ void rotateToLocal(const MeshLib::RotationMatrix& matR2local,
/// get a rotation matrix to the global coordinates
/// it computes R in x=R*x' where x is original coordinates and x' is local
/// coordinates
-void getRotationMatrixToGlobal(const unsigned element_dimension,
- const unsigned global_dim,
- const std::vector<MathLib::Point3d>& points,
- MeshLib::RotationMatrix& matR)
+MeshLib::RotationMatrix getRotationMatrixToGlobal(
+ const unsigned element_dimension,
+ const unsigned global_dim,
+ const std::vector<MathLib::Point3d>& points)
{
+ Eigen::Matrix3d matR;
// compute R in x=R*x' where x are original coordinates and x' are local
// coordinates
if (element_dimension == 1)
{
- auto const a = Eigen::Map<Eigen::Vector3d const>(points[0].getCoords());
- auto const b = Eigen::Map<Eigen::Vector3d const>(points[1].getCoords());
- Eigen::Vector3d xx = b - a;
- xx.normalize();
+ Eigen::Vector3d const xx =
+ (Eigen::Map<Eigen::Vector3d const>(points[1].getCoords()) -
+ Eigen::Map<Eigen::Vector3d const>(points[0].getCoords()))
+ .normalized();
if (global_dim == 2)
{
- GeoLib::compute2DRotationMatrixToX(xx, matR);
+ matR = GeoLib::compute2DRotationMatrixToX(xx);
}
else
{
- GeoLib::compute3DRotationMatrixToX(xx, matR);
+ matR = GeoLib::compute3DRotationMatrixToX(xx);
}
matR.transposeInPlace();
}
@@ -64,6 +65,7 @@ void getRotationMatrixToGlobal(const unsigned element_dimension,
// set a transposed matrix
matR.transposeInPlace();
}
+ return matR;
}
} // namespace detail
@@ -71,7 +73,7 @@ namespace MeshLib
{
ElementCoordinatesMappingLocal::ElementCoordinatesMappingLocal(
const Element& e, const unsigned global_dim)
- : _global_dim(global_dim), _matR2global(3, 3)
+ : _global_dim(global_dim), _matR2global(Eigen::Matrix3d::Identity())
{
assert(e.getDimension() <= global_dim);
_points.reserve(e.getNumberOfNodes());
@@ -82,14 +84,12 @@ ElementCoordinatesMappingLocal::ElementCoordinatesMappingLocal(
auto const element_dim = e.getDimension();
- if (global_dim == element_dim)
+ if (global_dim != element_dim)
{
- _matR2global.setIdentity();
- return;
+ _matR2global =
+ detail::getRotationMatrixToGlobal(element_dim, global_dim, _points);
+ detail::rotateToLocal(_matR2global.transpose(), _points);
}
-
- detail::getRotationMatrixToGlobal(element_dim, global_dim, _points, _matR2global);
- detail::rotateToLocal(_matR2global.transpose(), _points);
}
} // namespace MeshLib