diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h index cf5596402d65e07474a2d2628f06a7b696285612..8d8e91b58d42cef0a1e6e1592949134e8ffb9f33 100644 --- a/GeoLib/AnalyticalGeometry-impl.h +++ b/GeoLib/AnalyticalGeometry-impl.h @@ -70,59 +70,6 @@ std::pair getNewellPlane( return std::make_pair(plane_normal, d); } -template -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 -void compute3DRotationMatrixToX(Eigen::Vector3d const& v, T_MATRIX& rot_mat) -{ - double const eps = std::numeric_limits::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 void rotatePoints(Eigen::Matrix3d const& rot_mat, InputIterator pnts_begin, diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index cf481952d7ecfdf283dc9487673f3dcd8a52e904..61cafd20c91e413676dd6f0f872902c11666f291 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::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 3592d5c2f6d074c38fba7a502cf0663524bc99ab..46b4ebf987118310e7a868ce7fe36b01ffcf8c8d 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -78,20 +78,20 @@ std::pair getNewellPlane( const std::vector& 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 -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 -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 eb9fcb1e4afa44e738a7dc80ed6c56d3a38f1ff3..42aaaffdd77ebee4e3518b434dff92cb4868e5e1 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& points, - MeshLib::RotationMatrix& matR) +MeshLib::RotationMatrix getRotationMatrixToGlobal( + const unsigned element_dimension, + const unsigned global_dim, + const std::vector& 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(points[0].getCoords()); - auto const b = Eigen::Map(points[1].getCoords()); - Eigen::Vector3d xx = b - a; - xx.normalize(); + Eigen::Vector3d const xx = + (Eigen::Map(points[1].getCoords()) - + Eigen::Map(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