From 151e9d4400a74c728b470aaa25c2d6b1da746598 Mon Sep 17 00:00:00 2001 From: Norihiro Watanabe <norihiro.watanabe@ufz.de> Date: Mon, 27 Apr 2015 12:46:04 +0200 Subject: [PATCH] use a template parameter for a matrix type in computeRotationMatrixToXY() --- GeoLib/AnalyticalGeometry-impl.h | 58 ++++++++++++++++++++++++++++++++ GeoLib/AnalyticalGeometry.cpp | 57 ------------------------------- GeoLib/AnalyticalGeometry.h | 3 +- 3 files changed, 60 insertions(+), 58 deletions(-) diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h index e2ad36a5ec9..27adefe1b32 100644 --- a/GeoLib/AnalyticalGeometry-impl.h +++ b/GeoLib/AnalyticalGeometry-impl.h @@ -78,5 +78,63 @@ void compute3DRotationMatrixToX(MathLib::Vector3 const& v, } } +template <class T_MATRIX> +void computeRotationMatrixToXY(MathLib::Vector3 const& n, + T_MATRIX & rot_mat) +{ + // check if normal points already in the right direction + if (sqrt(n[0]*n[0]+n[1]*n[1]) == 0) { + 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) = 1.0; + rot_mat(1,2) = 0.0; + rot_mat(2,0) = 0.0; + rot_mat(2,1) = 0.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 h1(1 / n.getLength()); + + // 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; +} + } // end namespace GeoLib diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index bd9900cdd2e..67b5f1d6e0f 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -348,63 +348,6 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1, return std::abs(GeoLib::scalarTriple(ac, ad, ab)) / 6.0; } -void computeRotationMatrixToXY(MathLib::Vector3 const& n, - MathLib::DenseMatrix<double> & rot_mat) -{ - // check if normal points already in the right direction - if (sqrt(n[0]*n[0]+n[1]*n[1]) == 0) { - 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) = 1.0; - rot_mat(1,2) = 0.0; - rot_mat(2,0) = 0.0; - rot_mat(2,1) = 0.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 h1(1 / n.getLength()); - - // 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; -} - void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat) { // *** some frequently used terms *** diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index 3485bf32617..aea8be332bb 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -86,8 +86,9 @@ void compute3DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat); * @param plane_normal the (3d) vector that is rotated parallel to the \f$z\f$ axis * @param rot_mat 3x3 rotation matrix */ +template <class T_MATRIX> void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal, - MathLib::DenseMatrix<double> & rot_mat); + T_MATRIX & rot_mat); /** * Method computes the rotation matrix that rotates the given vector parallel to the \f$y\f$ axis. -- GitLab