Changeset 2050

Show
Ignore:
Timestamp:
08/28/08 14:13:05 (3 months ago)
Author:
bmsmith
Message:

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?.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • MOAB/trunk/tools/mcnpmit/main.cpp

    r1899 r2050  
    189189  std::vector<double> coords; 
    190190 
     191  double tal_sum     = 0.0, 
     192         err_sum     = 0.0, 
     193         tal_sum_sqr = 0.0, 
     194         err_sum_sqr = 0.0; 
     195 
    191196//  double davg = 0.0; 
    192197//  unsigned int    nmax = 0, nmin = 1000000000 ; 
    …… 
    221226    testvc[2] = transformed_pt[2]; 
    222227 
    223     // std::cout << n << " " << testvc << std::endl; 
    224  
    225228    // Find the leaf containing the point 
    226229    MBresult = kdtree.leaf_containing_point( root, transformed_pt, treeiter); 
    227     assert(MBresult == MB_SUCCESS); 
     230    if (MB_SUCCESS != MBresult) { 
     231      double x, y, z; 
     232      if (CARTESIAN == coord_sys) { 
     233        x = testvc[0]; 
     234        y = testvc[1]; 
     235        z = testvc[2]; 
     236      } 
     237      else if (CYLINDRICAL == coord_sys) { 
     238        x = testvc[0]*cos(2*M_PI*testvc[2]); 
     239        y = testvc[0]*sin(2*M_PI*testvc[2]); 
     240        z = testvc[1]; 
     241      }  
     242      else { 
     243        assert(MB_SUCCESS == MBresult); 
     244      } 
     245      std::cout << "No leaf found, MCNP coord xyz=" << x << " " << y << " " << z << std::endl; 
     246      cfd_iter++; 
     247      continue; 
     248    } 
    228249 
    229250    range.clear(); 
    …… 
    270291        found = true; 
    271292        elems_read++; 
     293         
     294        tal_sum     = tal_sum + taldata; 
     295        err_sum     = err_sum + errdata; 
     296        tal_sum_sqr = tal_sum_sqr + taldata*taldata; 
     297        err_sum_sqr = err_sum_sqr + errdata*errdata; 
    272298 
    273299          break; 
    …… 
    293319    std::cout << "Failure during query! " << elems_read << " elements interpolated." << std::endl << std::endl; 
    294320  } 
     321 
     322 
     323  double tal_std_dev = sqrt( (1.0/elems_read)*(tal_sum_sqr - (1.0/elems_read)*tal_sum*tal_sum) ); 
     324  double err_std_dev = sqrt( (1.0/elems_read)*(err_sum_sqr - (1.0/elems_read)*err_sum*err_sum) ); 
     325 
     326  std::cout << "Tally Mean:               " << tal_sum / elems_read << std::endl; 
     327  std::cout << "Tally Standard Deviation: " << tal_std_dev << std::endl; 
     328  std::cout << "Error Mean:               " << err_sum / elems_read << std::endl; 
     329  std::cout << "Error Standard Deviation: " << err_std_dev << std::endl; 
    295330 
    296331  interp_time = clock() - build_time; 

  • MOAB/trunk/tools/mcnpmit/mcnpmit.cpp

    r1899 r2050  
    422422      switch( csys ) { 
    423423        case CARTESIAN : 
    424           r[0] = q[0]; r[1] = q[1]; r[2] = q[2]; 
     424          r[0] = q[0]; r[1] = q[1]; r[2] = q[2];  // x, y, z 
    425425        break; 
    426426        case CYLINDRICAL : 
    427           r[0] = sqrt( q[0]*q[0] + q[1]*q[1] ); 
    428           r[1] = q[2]; 
    429           r[2] = c2pi * ( atan2( q[1], q[0] ) ); 
     427          r[0] = sqrt( q[0]*q[0] + q[1]*q[1] );   // r 
     428          r[1] = q[2];                            // z 
     429          r[2] = c2pi * ( atan2( q[1], q[0] ) );  // theta (in rotations) 
    430430        break; 
    431431        case SPHERICAL :