MOAB: Mesh Oriented datABase  (version 5.4.1)
main.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <cmath>
00004 #include <ctime>
00005 #include <cstdlib>
00006 #include <cassert>
00007 
00008 #include "mcnpmit.hpp"
00009 
00010 #include "moab/CartVect.hpp"
00011 #include "moab/Core.hpp"
00012 #include "MBTagConventions.hpp"
00013 #include "moab/AdaptiveKDTree.hpp"
00014 #include "moab/GeomUtil.hpp"
00015 #include "moab/FileOptions.hpp"
00016 #include "../tools/mbcoupler/ElemUtil.hpp"
00017 
00018 #define MBI mb_instance()
00019 
00020 McnpData* mc_instance();
00021 moab::Interface* mb_instance();
00022 MCNPError read_files( int, char** );
00023 MCNPError next_double( std::string, double&, int& );
00024 MCNPError next_int( std::string, int&, int& );
00025 
00026 moab::Tag coord_tag, rotation_tag, cfd_heating_tag, cfd_error_tag;
00027 
00028 std::string h5m_filename;
00029 std::string CAD_filename;
00030 std::string output_filename;
00031 
00032 bool skip_build = false;
00033 bool read_qnv   = false;
00034 
00035 clock_t start_time, load_time, build_time, interp_time;
00036 
00037 int main( int argc, char** argv )
00038 {
00039     MCNPError result;
00040 
00041     start_time = clock();
00042 
00043     // Read in file names from command line
00044     result = read_files( argc, argv );
00045     if( result == MCNP_FAILURE ) return 1;
00046 
00047     result = MCNP->initialize_tags();
00048 
00049     // Parse the MCNP input file
00050     if( !skip_build )
00051     {
00052 
00053         result = MCNP->read_mcnpfile( skip_build );
00054         if( result == MCNP_FAILURE )
00055         {
00056             std::cout << "Failure reading MCNP file!" << std::endl;
00057             return 1;
00058         }
00059     }
00060 
00061     load_time = clock() - start_time;
00062 
00063     // Make the KD-Tree
00064     moab::ErrorCode MBresult;
00065     moab::AdaptiveKDTree kdtree( MBI );
00066     moab::EntityHandle root;
00067 
00068     MBI->tag_get_handle( "CoordTag", 1, moab::MB_TYPE_INTEGER, coord_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
00069     MBI->tag_get_handle( "RotationTag", 16, moab::MB_TYPE_DOUBLE, rotation_tag,
00070                          moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
00071 
00072     if( skip_build )
00073     {
00074         MBresult = MBI->load_mesh( h5m_filename.c_str() );
00075 
00076         if( moab::MB_SUCCESS == MBresult )
00077         {
00078             std::cout << std::endl << "Read in mesh from h5m file." << std::endl << std::endl;
00079             std::cout << "Querying mesh file..." << std::endl;
00080         }
00081         else
00082         {
00083             std::cout << "Failure reading h5m file!" << std::endl;
00084             std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
00085             std::string message;
00086             if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
00087                 std::cerr << "Error message: " << message << std::endl;
00088             return 1;
00089         }
00090 
00091         moab::Range tmprange;
00092         kdtree.find_all_trees( tmprange );
00093         root = tmprange[0];
00094     }
00095     else
00096     {
00097         std::cout << "Building KD-Tree..." << std::endl;
00098         moab::FileOptions opts( "CANDIDATE_PLANE_SET=SUBDIVISION" );
00099         MBresult = kdtree.build_tree( MCNP->elem_handles, &root, &opts );
00100         if( MBresult == moab::MB_SUCCESS )
00101         {
00102 
00103             MBI->tag_set_data( coord_tag, &root, 1, &( MCNP->coord_system ) );
00104             MBI->tag_set_data( rotation_tag, &root, 1, &( MCNP->rotation_matrix ) );
00105 
00106             std::cout << "KD-Tree has been built successfully!" << std::endl << std::endl;
00107             MBresult = MBI->write_mesh( ( MCNP->MCNP_filename + ".h5m" ).c_str() );
00108 
00109             std::cout << "Querying mesh file..." << std::endl;
00110         }
00111         else
00112         {
00113             std::cout << "Error building KD-Tree!" << std::endl << std::endl;
00114             std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
00115             std::string message;
00116             if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
00117                 std::cerr << "Error message: " << message << std::endl;
00118             return 1;
00119         }
00120     }
00121 
00122     int coord_sys;
00123     double rmatrix[16];
00124 
00125     MBresult = MBI->tag_get_data( coord_tag, &root, 1, &coord_sys );MB_CHK_ERR( MBresult );
00126     MBresult = MBI->tag_get_data( rotation_tag, &root, 1, &rmatrix );MB_CHK_ERR( MBresult );
00127 
00128     build_time = clock() - load_time;
00129 
00130     // Read the CAD mesh data and query the tree
00131     std::ifstream cadfile;
00132     std::ofstream outfile;
00133 
00134     outfile.open( output_filename.c_str() );
00135 
00136     int num_pts;
00137     int n;
00138     long int nl = 0;
00139     char* ctmp;
00140     int elems_read = 0;
00141     int p          = 0;
00142     char line[10000];
00143 
00144     // Used only when reading a mesh file to get vertex info
00145     double* cfd_coords = NULL;
00146     moab::Range::iterator cfd_iter;
00147     moab::EntityHandle meshset;
00148 
00149     if( read_qnv )
00150     {
00151         cadfile.open( CAD_filename.c_str() );
00152         cadfile.getline( line, 10000 );
00153         cadfile.getline( line, 10000 );
00154         result = next_int( line, num_pts, p );
00155     }
00156     else
00157     {
00158 
00159         meshset  = 0;
00160         MBresult = MBI->load_file( CAD_filename.c_str(), &meshset );MB_CHK_ERR( MBresult );
00161         assert( 0 != meshset );
00162 
00163         moab::Range cfd_verts;
00164         MBresult = MBI->get_entities_by_type( meshset, moab::MBVERTEX, cfd_verts, true );MB_CHK_ERR( MBresult );
00165         num_pts = cfd_verts.size();
00166 
00167         cfd_coords = new double[3 * num_pts];
00168         MBresult   = MBI->get_coords( cfd_verts, cfd_coords );MB_CHK_ERR( MBresult );
00169 
00170         cfd_iter = cfd_verts.begin();
00171         MBresult = MBI->tag_get_handle( "heating_tag", 1, moab::MB_TYPE_DOUBLE, cfd_heating_tag,
00172                                         moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );
00173         MBresult = MBI->tag_get_handle( "error_tag", 1, moab::MB_TYPE_DOUBLE, cfd_error_tag,
00174                                         moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );
00175 
00176         std::cout << std::endl << "Read in mesh with query points." << std::endl << std::endl;
00177     }
00178 
00179     double testpt[3];
00180     double transformed_pt[3];
00181     double taldata;
00182     double errdata;
00183 
00184     moab::CartVect testvc;
00185 
00186     bool found = false;
00187 
00188     // MBRange verts;
00189     std::vector< moab::EntityHandle > verts;
00190     moab::Range range;
00191     moab::CartVect box_max, box_min;
00192 
00193     moab::CartVect hexverts[8];
00194     moab::CartVect tmp_cartvect;
00195     std::vector< double > coords;
00196 
00197     double tal_sum = 0.0, err_sum = 0.0, tal_sum_sqr = 0.0, err_sum_sqr = 0.0;
00198 
00199     //  double davg = 0.0;
00200     //  unsigned int    nmax = 0, nmin = 1000000000 ;
00201 
00202     for( unsigned int i = 0; i < (unsigned int)num_pts; i++ )
00203     {
00204 
00205         // if (i%status_freq == 0)
00206         //  std::cerr << "Completed " << i/status_freq << "%" << std::endl;
00207 
00208         // Grab the coordinates to query
00209         if( read_qnv )
00210         {
00211             cadfile.getline( line, 10000 );
00212 
00213             nl        = std::strtol( line, &ctmp, 10 );
00214             n         = (unsigned int)nl;
00215             testpt[0] = std::strtod( ctmp + 1, &ctmp );
00216             testpt[1] = std::strtod( ctmp + 1, &ctmp );
00217             testpt[2] = std::strtod( ctmp + 1, NULL );
00218         }
00219         else
00220         {
00221             testpt[0] = cfd_coords[3 * i];
00222             testpt[1] = cfd_coords[3 * i + 1];
00223             testpt[2] = cfd_coords[3 * i + 2];
00224             n         = i + 1;
00225         }
00226 
00227         result = MCNP->transform_point( testpt, transformed_pt, coord_sys, rmatrix );
00228 
00229         testvc[0] = transformed_pt[0];
00230         testvc[1] = transformed_pt[1];
00231         testvc[2] = transformed_pt[2];
00232 
00233         // Find the leaf containing the point
00234         moab::EntityHandle tree_node;
00235         MBresult = kdtree.point_search( transformed_pt, tree_node );
00236         if( moab::MB_SUCCESS != MBresult )
00237         {
00238             double x = 0.0, y = 0.0, z = 0.0;
00239             if( CARTESIAN == coord_sys )
00240             {
00241                 x = testvc[0];
00242                 y = testvc[1];
00243                 z = testvc[2];
00244             }
00245             else if( CYLINDRICAL == coord_sys )
00246             {
00247                 x = testvc[0] * cos( 2 * M_PI * testvc[2] );
00248                 y = testvc[0] * sin( 2 * M_PI * testvc[2] );
00249                 z = testvc[1];
00250             }
00251             else
00252             {
00253                 std::cout << "MOAB WARNING: Unhandled error code during point search in KdTree, "
00254                              "ErrorCode = "
00255                           << MBresult << " and Coord xyz=" << x << " " << y << " " << z << std::endl;
00256             }
00257             std::cout << "No leaf found, MCNP coord xyz=" << x << " " << y << " " << z << std::endl;
00258             ++cfd_iter;
00259             continue;
00260         }
00261 
00262         range.clear();
00263         MBresult = MBI->get_entities_by_type( tree_node, moab::MBHEX, range );MB_CHK_ERR( MBresult );
00264 
00265         // davg += (double) range.size();
00266         // if (range.size() > nmax) nmax = range.size();
00267         // if (range.size() < nmin) nmin = range.size();
00268 
00269         for( moab::Range::iterator rit = range.begin(); rit != range.end(); ++rit )
00270         {
00271             verts.clear();
00272             const moab::EntityHandle* connect;
00273             int num_connect;
00274             MBresult = MBI->get_connectivity( *rit, connect, num_connect, true );MB_CHK_ERR( MBresult );
00275 
00276             coords.resize( 3 * num_connect );
00277             MBresult = MBI->get_coords( connect, num_connect, &coords[0] );MB_CHK_ERR( MBresult );
00278 
00279             for( unsigned int j = 0; j < (unsigned int)num_connect; j++ )
00280             {
00281                 hexverts[j][0] = coords[3 * j];
00282                 hexverts[j][1] = coords[3 * j + 1];
00283                 hexverts[j][2] = coords[3 * j + 2];
00284             }
00285 
00286             if( moab::ElemUtil::point_in_trilinear_hex( hexverts, testvc, 1.e-6 ) )
00287             {
00288                 MBresult = MBI->tag_get_data( MCNP->tally_tag, &( *rit ), 1, &taldata );MB_CHK_ERR( MBresult );
00289                 MBresult = MBI->tag_get_data( MCNP->relerr_tag, &( *rit ), 1, &errdata );MB_CHK_ERR( MBresult );
00290 
00291                 outfile << n << "," << testpt[0] << "," << testpt[1] << "," << testpt[2] << "," << taldata << ","
00292                         << errdata << std::endl;
00293 
00294                 if( !read_qnv )
00295                 {
00296                     MBresult = MBI->tag_set_data( cfd_heating_tag, &( *cfd_iter ), 1, &taldata );MB_CHK_ERR( MBresult );
00297                     MBresult = MBI->tag_set_data( cfd_error_tag, &( *cfd_iter ), 1, &errdata );MB_CHK_ERR( MBresult );
00298                 }
00299 
00300                 found = true;
00301                 elems_read++;
00302 
00303                 tal_sum     = tal_sum + taldata;
00304                 err_sum     = err_sum + errdata;
00305                 tal_sum_sqr = tal_sum_sqr + taldata * taldata;
00306                 err_sum_sqr = err_sum_sqr + errdata * errdata;
00307 
00308                 break;
00309             }
00310         }
00311 
00312         if( !read_qnv ) ++cfd_iter;
00313 
00314         if( !found )
00315         {
00316             std::cout << n << " " << testvc << std::endl;
00317         }
00318         found = false;
00319     }
00320 
00321     cadfile.close();
00322     outfile.close();
00323 
00324     if( result == MCNP_SUCCESS )
00325     {
00326         std::cout << "Success! " << elems_read << " elements interpolated." << std::endl << std::endl;
00327     }
00328     else
00329     {
00330         std::cout << "Failure during query! " << elems_read << " elements interpolated." << std::endl << std::endl;
00331     }
00332 
00333     double tal_std_dev = sqrt( ( 1.0 / elems_read ) * ( tal_sum_sqr - ( 1.0 / elems_read ) * tal_sum * tal_sum ) );
00334     double err_std_dev = sqrt( ( 1.0 / elems_read ) * ( err_sum_sqr - ( 1.0 / elems_read ) * err_sum * err_sum ) );
00335 
00336     std::cout << "Tally Mean:               " << tal_sum / elems_read << std::endl;
00337     std::cout << "Tally Standard Deviation: " << tal_std_dev << std::endl;
00338     std::cout << "Error Mean:               " << err_sum / elems_read << std::endl;
00339     std::cout << "Error Standard Deviation: " << err_std_dev << std::endl;
00340 
00341     interp_time = clock() - build_time;
00342 
00343     if( !read_qnv )
00344     {
00345         std::string out_mesh_fname = output_filename;
00346         MBresult                   = MBI->write_mesh( ( out_mesh_fname + ".h5m" ).c_str(), &meshset, 1 );
00347         // MBresult = MBI->write_file( (cfd_mesh_fname + ".vtk").c_str(), "vtk", NULL, &meshset, 1,
00348         // &cfd_heating_tag, 1);
00349     }
00350 
00351     std::cout << "Time to read in file:     " << (double)load_time / CLOCKS_PER_SEC << std::endl;
00352     std::cout << "Time to build kd-tree:    " << (double)build_time / CLOCKS_PER_SEC << std::endl;
00353     std::cout << "Time to interpolate data: " << (double)interp_time / CLOCKS_PER_SEC << std::endl;
00354 
00355     return 0;
00356 }
00357 
00358 MCNPError read_files( int argc, char** argv )
00359 {
00360     MCNPError result = MCNP_FAILURE;
00361 
00362     // Check to see if appropriate command lines specified
00363     if( argc < 3 )
00364     {
00365         std::cout << "Source and Target mesh filenames NOT specified!";
00366         std::cout << std::endl;
00367         return MCNP_FAILURE;
00368     }
00369 
00370     // Set the MCNP or H5M filename
00371     std::string str;
00372     str = argv[1];
00373 
00374     unsigned int itmp = str.find( ".h5m" );
00375     if( ( itmp > 0 ) && ( itmp < str.length() ) )
00376     {
00377         skip_build   = true;
00378         h5m_filename = str;
00379     }
00380     else
00381     {
00382         result = MCNP->set_filename( str );
00383     }
00384 
00385     // Set the CAD filename
00386     str          = argv[2];
00387     CAD_filename = str;
00388 
00389     itmp = str.find( ".qnv" );
00390     if( ( itmp > 0 ) && ( itmp < str.length() ) ) read_qnv = true;
00391 
00392     // Set the output filename
00393     str             = argv[3];
00394     output_filename = str;
00395 
00396     return result;
00397 }
00398 
00399 MCNPError next_double( std::string s, double& d, int& p )
00400 {
00401 
00402     unsigned int slen = s.length();
00403     unsigned int j;
00404 
00405     for( unsigned int i = p; i < slen; i++ )
00406     {
00407         if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
00408         {
00409 
00410             j = s.find( ",", i );
00411             if( j > slen ) j = slen;
00412 
00413             d = std::atof( s.substr( i, j - i ).c_str() );
00414             p = j + 1;
00415 
00416             return MCNP_SUCCESS;
00417         }
00418     }
00419 
00420     return DONE;
00421 }
00422 
00423 MCNPError next_int( std::string s, int& k, int& p )
00424 {
00425 
00426     unsigned int slen = s.length();
00427     unsigned int j;
00428 
00429     for( unsigned int i = p; i < slen; i++ )
00430     {
00431         if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
00432         {
00433 
00434             j = s.find( ",", i );
00435             if( j > slen ) j = slen;
00436 
00437             k = std::atoi( s.substr( i, j - i ).c_str() );
00438             p = j + 1;
00439 
00440             return MCNP_SUCCESS;
00441         }
00442     }
00443 
00444     return DONE;
00445 }
00446 
00447 McnpData* mc_instance()
00448 {
00449     static McnpData inst;
00450     return &inst;
00451 }
00452 
00453 moab::Interface* mb_instance()
00454 {
00455     static moab::Core inst;
00456     return &inst;
00457 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines