diff --git a/src/OrientedBox.cpp b/src/OrientedBox.cpp index a4d5f358d872c86946fb02882ac6a887a00322df..d90c24b8fa5451a4036aa1eb200b3bd7d6721e80 100644 --- a/src/OrientedBox.cpp +++ b/src/OrientedBox.cpp @@ -563,9 +563,8 @@ bool OrientedBox::intersect_ray( const CartVect& ray_origin, } // get transpose of axes - Matrix3 B = Matrix::transpose(axes); + Matrix3 B = axes.transpose(); - // transform ray to box coordintae system CartVect par_pos = B * (ray_origin - center); CartVect par_dir = B * ray_direction; diff --git a/src/moab/Matrix3.hpp b/src/moab/Matrix3.hpp index 37f172fb9d9318eb82431fcb9535e4d1fd7a81b0..fb212abc9af583487a14486bcc21bf45927d0d17 100644 --- a/src/moab/Matrix3.hpp +++ b/src/moab/Matrix3.hpp @@ -55,17 +55,17 @@ namespace moab { namespace Matrix{ template< typename Matrix> - Matrix inverse( const Matrix & d, const double i){ + Matrix inverse( const Matrix & d, const double det){ Matrix m( d); - m( 0) = i * (d(4) * d(8) - d(5) * d(7)); - m( 1) = i * (d(2) * d(7) - d(8) * d(1)); - m( 2) = i * (d(1) * d(5) - d(4) * d(2)); - m( 3) = i * (d(5) * d(6) - d(8) * d(3)); - m( 4) = i * (d(0) * d(8) - d(6) * d(2)); - m( 5) = i * (d(2) * d(3) - d(5) * d(0)); - m( 6) = i * (d(3) * d(7) - d(6) * d(4)); - m( 7) = i * (d(1) * d(6) - d(7) * d(0)); - m( 8) = i * (d(0) * d(4) - d(3) * d(1)); + m(0) = det * (d(4) * d(8) - d(5) * d(7)); + m(1) = det * (d(2) * d(7) - d(8) * d(1)); + m(2) = det * (d(1) * d(5) - d(4) * d(2)); + m(3) = det * (d(5) * d(6) - d(8) * d(3)); + m(4) = det * (d(0) * d(8) - d(6) * d(2)); + m(5) = det * (d(2) * d(3) - d(5) * d(0)); + m(6) = det * (d(3) * d(7) - d(6) * d(4)); + m(7) = det * (d(1) * d(6) - d(7) * d(0)); + m(8) = det * (d(0) * d(4) - d(3) * d(1)); return m; } @@ -78,12 +78,7 @@ namespace Matrix{ det = d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8; return d(0) > 0 && subdet8 > 0 && det > 0; } - template< typename Matrix> - inline Matrix transpose( const Matrix & d){ - return Matrix( d(0), d(3), d(6), - d(1), d(4), d(7), - d(2), d(5), d(8) ); - } + template< typename Matrix> inline Matrix mmult3( const Matrix& a, const Matrix& b ) { return Matrix( a(0,0) * b(0,0) + a(0,1) * b(1,0) + a(0,2) * b(2,0), @@ -138,150 +133,6 @@ namespace Matrix{ return res; } - // moved from EigenDecomp.hpp - - // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn - // real symmetric matrix. Square nxn matrix a; size of matrix in n; - // output eigenvalues in w; and output eigenvectors in v. Resulting - // eigenvalues/vectors are sorted in decreasing order; eigenvectors are - // normalized. - // TODO: Remove this - #define VTK_ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ - a[k][l]=h+s*(g-h*tau) - - //TODO: Refactor this method into subroutines - //use a namespace { } with no name to - //contain subroutines so that the compiler - //automatically inlines them. - - template< typename Matrix, typename Vector> - ErrorCode EigenDecomp2( const Matrix & _a, - double w[3], - Vector o[3] ) { - Vector v[3]; - const int MAX_ROTATIONS = 20; - const double one_ninth = 1./9; - int i, j, k, iq, ip, numPos; - double tresh, theta, tau, t, sm, s, h, g, c, tmp; - double b[3], z[3]; - Matrix a( _a); - - // initialize - for (ip=0; ip<3; ip++) { - for (iq=0; iq<3; iq++){ - v[ip][iq] = 0.0; - } - v[ip][ip] = 1.0; - } - for (ip=0; ip<3; ip++) { - b[ip] = w[ip] = a[ip][ip]; - z[ip] = 0.0; - } - - // begin rotation sequence - for (i=0; i 3 && (fabs(w[ip])+g) == fabs(w[ip]) - && (fabs(w[iq])+g) == fabs(w[iq])) { - a[ip][iq] = 0.0; - } - else if ( fabs(a[ip][iq]) > tresh) { - h = w[iq] - w[ip]; - if ( (fabs(h)+g) == fabs(h)){ t = (a[ip][iq]) / h; } - else { - theta = 0.5*h / (a[ip][iq]); - t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta)); - if (theta < 0.0) { t = -t;} - } - c = 1.0 / sqrt(1+t*t); - s = t*c; - tau = s/(1.0+c); - h = t*a[ip][iq]; - z[ip] -= h; - z[iq] += h; - w[ip] -= h; - w[iq] += h; - a[ip][iq]=0.0; - // ip already shifted left by 1 unit - for (j = 0;j <= ip-1;j++) { VTK_ROTATE(a,j,ip,j,iq); } - // ip and iq already shifted left by 1 unit - for (j = ip+1;j <= iq-1;j++) { VTK_ROTATE(a,ip,j,j,iq); } - // iq already shifted left by 1 unit - for (j=iq+1; j<3; j++) { VTK_ROTATE(a,ip,j,iq,j); } - for (j=0; j<3; j++) { VTK_ROTATE(v,j,ip,j,iq); } - } - } - } - - for (ip=0; ip<3; ip++) { - b[ip] += z[ip]; - w[ip] = b[ip]; - z[ip] = 0.0; - } - } - - //// this is NEVER called - if ( i >= MAX_ROTATIONS ) { - std::cerr << "Matrix3D: Error extracting eigenfunctions" << std::endl; - return MB_FAILURE; - } - - // sort eigenfunctions these changes do not affect accuracy - for (j=0; j<2; j++){ // boundary incorrect - k = j; - tmp = w[k]; - for (i=j+1; i<3; i++){ // boundary incorrect, shifted already - if (w[i] >= tmp){ // why exchage if same? - k = i; - tmp = w[k]; - } - } - if (k != j){ - w[k] = w[j]; - w[j] = tmp; - for (i=0; i<3; i++){ - tmp = v[i][j]; - v[i][j] = v[i][k]; - v[i][k] = tmp; - } - } - } - // insure eigenvector consistency (i.e., Jacobi can compute vectors that - // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can - // reek havoc in hyperstreamline/other stuff. We will select the most - // positive eigenvector. - int ceil_half_n = (3 >> 1) + (3 & 1); - for (j=0; j<3; j++) { - for (numPos=0, i=0; i<3; i++) { - if ( v[i][j] >= 0.0 ) { numPos++; } - } - // if ( numPos < ceil(double(n)/double(2.0)) ) - if ( numPos < ceil_half_n) { - for(i=0; i<3; i++) { v[i][j] *= -1.0; } - } - } - - //transpose the vector array for output - for(i=0; i<3; i++) - { - for(j=0; j<3; j++) - { o[i][j]=v[j][i]; } - } - - return MB_SUCCESS; - } } //namespace Matrix class Matrix3 {