Merged in iulian07/matrix3_rows (pull request #266)

Matrix3 constructor with 3 rows

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;

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<MAX_ROTATIONS; i++){
-	    sm = 0.0;
-	    for (ip=0; ip<2; ip++){
-	      for (iq=ip+1; iq<3; iq++){ sm += fabs(a[ip][iq]); }
-	    }
-
-	    if ( sm == 0.0 ){ break; }
-	    // first 3 sweeps
-	    tresh = (i < 3)? 0.2*sm*one_ninth : 0.0;
-	    for (ip=0; ip<2; ip++) {
-	      for (iq=ip+1; iq<3; iq++) {
-	        g = 100.0*fabs(a[ip][iq]);
-
-	        // after 4 sweeps
-	        if ( 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  {

tools/dagmc/DagMC.cpp

@@ -65,7 +65,7 @@ namespace moab {
 const std::map<std::string, std::string> DagMC::no_synonyms;
 
 // DagMC Constructor
-DagMC::DagMC(Interface *mb_impl, double overlap_tolerance, double numerical_precision) {
+DagMC::DagMC(Interface *mb_impl, double overlap_tolerance, double p_numerical_precision) {
   moab_instance_created = false;
   // if we arent handed a moab instance create one
   if (NULL == mb_impl) {
@@ -82,7 +82,7 @@ DagMC::DagMC(Interface *mb_impl, double overlap_tolerance, double numerical_prec
   // This is the correct place to uniquely define default values for the dagmc settings
   overlapThickness = overlap_tolerance; // must be nonnegative
   defaultFacetingTolerance = .001;
-  numericalPrecision = numerical_precision;
+  numericalPrecision = p_numerical_precision;
 
   memset( implComplName, 0, NAME_TAG_SIZE );
   strcpy( implComplName , "impl_complement" );

tools/mcnpmit/main.cpp

@@ -23,7 +23,6 @@ MCNPError read_files(int, char**);
 MCNPError next_double(std::string, double&, int&);
 MCNPError next_int(std::string, int&, int&);
 
-MCNPError result;
 moab::Tag     coord_tag, rotation_tag, cfd_heating_tag, cfd_error_tag;
 
 std::string h5m_filename;
@@ -36,6 +35,7 @@ bool read_qnv   = false;
 clock_t start_time, load_time, build_time, interp_time;
 
 int main(int argc, char **argv) {
+  MCNPError result;
 
   start_time = clock();
 
@@ -340,6 +340,7 @@ int main(int argc, char **argv) {
 
 
 MCNPError read_files(int argc, char **argv) {
+  MCNPError result;
 
   // Check to see if appropriate command lines specified
   if (argc < 3) {