[MOAB-dev] r2050 - MOAB/trunk/tools/mcnpmit

[bmsmith at mcs.anl.gov]

Author: bmsmith
Date: 2008-08-28 14:13:05 -0500 (Thu, 28 Aug 2008)
New Revision: 2050

Modified:
   MOAB/trunk/tools/mcnpmit/main.cpp
   MOAB/trunk/tools/mcnpmit/mcnpmit.cpp
Log:
Sometimes a cfd point cannot be found in the tree. Instead of an assert failing, I've added an error message.
Added calculation of error and standard deviation, possibly useful for setting plot range in VisIt.


Modified: MOAB/trunk/tools/mcnpmit/main.cpp
===================================================================
--- MOAB/trunk/tools/mcnpmit/main.cpp	2008-08-26 03:29:52 UTC (rev 2049)
+++ MOAB/trunk/tools/mcnpmit/main.cpp	2008-08-28 19:13:05 UTC (rev 2050)
@@ -188,6 +188,11 @@
   MBCartVect tmp_cartvect;
   std::vector<double> coords;
 
+  double tal_sum     = 0.0,
+         err_sum     = 0.0,
+         tal_sum_sqr = 0.0,
+         err_sum_sqr = 0.0;
+
 //  double davg = 0.0;
 //  unsigned int    nmax = 0, nmin = 1000000000 ;
 
@@ -220,11 +225,27 @@
     testvc[1] = transformed_pt[1]; 
     testvc[2] = transformed_pt[2];
 
-    // std::cout << n << " " << testvc << std::endl;
-
     // Find the leaf containing the point
     MBresult = kdtree.leaf_containing_point( root, transformed_pt, treeiter);
-    assert(MBresult == MB_SUCCESS);
+    if (MB_SUCCESS != MBresult) {
+      double x, y, z;
+      if (CARTESIAN == coord_sys) {
+        x = testvc[0];
+        y = testvc[1];
+	z = testvc[2];
+      }
+      else if (CYLINDRICAL == coord_sys) {
+        x = testvc[0]*cos(2*M_PI*testvc[2]);
+        y = testvc[0]*sin(2*M_PI*testvc[2]);
+        z = testvc[1];
+      } 
+      else {
+        assert(MB_SUCCESS == MBresult);
+      }
+      std::cout << "No leaf found, MCNP coord xyz=" << x << " " << y << " " << z << std::endl;
+      cfd_iter++;
+      continue;
+    }
 
     range.clear();
     MBresult = MBI -> get_entities_by_type( treeiter.handle(), MBHEX, range );
@@ -269,6 +290,11 @@
 
         found = true;
         elems_read++;
+	
+	tal_sum     = tal_sum + taldata;
+	err_sum     = err_sum + errdata;
+	tal_sum_sqr = tal_sum_sqr + taldata*taldata;
+	err_sum_sqr = err_sum_sqr + errdata*errdata;
 
 	  break;
       }
@@ -293,6 +319,15 @@
     std::cout << "Failure during query! " << elems_read << " elements interpolated." << std::endl << std::endl;
   }
 
+
+  double tal_std_dev = sqrt( (1.0/elems_read)*(tal_sum_sqr - (1.0/elems_read)*tal_sum*tal_sum) );
+  double err_std_dev = sqrt( (1.0/elems_read)*(err_sum_sqr - (1.0/elems_read)*err_sum*err_sum) );
+
+  std::cout << "Tally Mean:               " << tal_sum / elems_read << std::endl;
+  std::cout << "Tally Standard Deviation: " << tal_std_dev << std::endl;
+  std::cout << "Error Mean:               " << err_sum / elems_read << std::endl;
+  std::cout << "Error Standard Deviation: " << err_std_dev << std::endl;
+
   interp_time = clock() - build_time;
 
   if (!read_qnv) {

Modified: MOAB/trunk/tools/mcnpmit/mcnpmit.cpp
===================================================================
--- MOAB/trunk/tools/mcnpmit/mcnpmit.cpp	2008-08-26 03:29:52 UTC (rev 2049)
+++ MOAB/trunk/tools/mcnpmit/mcnpmit.cpp	2008-08-28 19:13:05 UTC (rev 2050)
@@ -421,12 +421,12 @@
       // Transform coordinate system
       switch( csys ) {
         case CARTESIAN :
-          r[0] = q[0]; r[1] = q[1]; r[2] = q[2];
+          r[0] = q[0]; r[1] = q[1]; r[2] = q[2];  // x, y, z
         break;
         case CYLINDRICAL :
-          r[0] = sqrt( q[0]*q[0] + q[1]*q[1] );
-          r[1] = q[2];
-          r[2] = c2pi * ( atan2( q[1], q[0] ) );
+          r[0] = sqrt( q[0]*q[0] + q[1]*q[1] );   // r
+          r[1] = q[2];                            // z
+          r[2] = c2pi * ( atan2( q[1], q[0] ) );  // theta (in rotations)
         break;
         case SPHERICAL :
           return MCNP_FAILURE;