Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
mcnpmit.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <cstdlib>
00004 #include "mcnpmit.hpp"
00005 #include "moab/CartVect.hpp"
00006 #include <cmath>
00007 
00008 moab::Interface* mb_instance();
00009 
00010 // Parameters
00011 // const double pi   = 3.141592653589793;
00012 const double c2pi = 0.1591549430918954;
00013 // const double cpi  = 0.3183098861837907;
00014 
00015 MCNPError next_number( std::string, double&, int& );
00016 int how_many_numbers( std::string );
00017 MCNPError read_numbers( std::string, int, std::vector< double >& );
00018 
00019 // Constructor
00020 McnpData::McnpData()
00021 {
00022 
00023     // Default value for coordinate system
00024     coord_system = 0;
00025 
00026     // Default rotation matrix is identity matrix
00027     for( int i = 0; i < 4; i++ )
00028     {
00029         for( int j = 0; j < 4; j++ )
00030         {
00031             if( i == j )
00032                 rotation_matrix[4 * i + j] = 1;
00033             else
00034                 rotation_matrix[4 * i + j] = 0;
00035         }
00036     }
00037 }
00038 
00039 // Destructor
00040 McnpData::~McnpData()
00041 {
00042 
00043     // Vertices and elements
00044     MCNP_vertices.clear();
00045 }
00046 
00047 // Setting and retrieving coordinate sysem
00048 MCNPError McnpData::set_coord_system( int k )
00049 {
00050     coord_system = k;
00051     return MCNP_SUCCESS;
00052 }
00053 int McnpData::get_coord_system()
00054 {
00055     return coord_system;
00056 }
00057 
00058 // Setting and retrieving roation matrix
00059 MCNPError McnpData::set_rotation_matrix( double r[16] )
00060 {
00061     for( int i = 0; i < 16; i++ )
00062     {
00063         rotation_matrix[i] = r[i];
00064     }
00065     return MCNP_SUCCESS;
00066 }
00067 double* McnpData::get_rotation_matrix()
00068 {
00069     return rotation_matrix;
00070 }
00071 
00072 // Set the filename
00073 MCNPError McnpData::set_filename( std::string fname )
00074 {
00075     MCNP_filename = fname;
00076     return MCNP_SUCCESS;
00077 }
00078 std::string McnpData::get_filename()
00079 {
00080     return MCNP_filename;
00081 }
00082 
00083 // Reading the MCNP file
00084 MCNPError McnpData::read_mcnpfile( bool skip_mesh )
00085 {
00086 
00087     MCNPError result;
00088     moab::ErrorCode MBresult;
00089     moab::CartVect tvect;
00090 
00091     std::vector< double > xvec[3];
00092 
00093     // Open the file
00094     std::ifstream mcnpfile;
00095     mcnpfile.open( MCNP_filename.c_str() );
00096     if( !mcnpfile )
00097     {
00098         std::cout << "Unable to open MCNP data file." << std::endl;
00099         return MCNP_FAILURE;
00100     }
00101     std::cout << std::endl;
00102     std::cout << "Reading MCNP input file..." << std::endl;
00103 
00104     // Prepare for file reading ...
00105     char line[10000];
00106     int mode = 0;  // Set the file reading mode to read proper data
00107     int nv[3];
00108 
00109     // Read in the file ...
00110     while( !mcnpfile.eof() )
00111     {
00112 
00113         mcnpfile.getline( line, 10000 );
00114         // std::cout << line << std::endl;
00115 
00116         switch( mode )
00117         {
00118             case 0:  // First line is a title
00119                 mode++;
00120                 break;
00121             case 1:  // Coordinate system
00122                 mode++;
00123                 result = read_coord_system( line );
00124                 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00125                 break;
00126             case 2:  // Rotation matrix
00127                 mode++;
00128                 for( int i = 0; i < 4; i++ )
00129                 {
00130                     mcnpfile.getline( line, 10000 );
00131                     result = read_rotation_matrix( line, i );
00132                     if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00133                 }
00134                 if( skip_mesh ) return MCNP_SUCCESS;
00135                 break;
00136             case 3:  // Read in vertices and build elements
00137                 mode++;
00138 
00139                 for( int i = 0; i < 3; i++ )
00140                 {
00141                     // How many points in the x[i]-direction
00142                     nv[i] = how_many_numbers( line );
00143                     if( nv[i] <= 0 ) return MCNP_FAILURE;
00144 
00145                     // Get space and read in these points
00146                     result = read_numbers( line, nv[i], xvec[i] );
00147                     if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00148 
00149                     // Update to the next line
00150                     mcnpfile.getline( line, 10000 );
00151                 }
00152 
00153                 // Make the elements and vertices
00154                 result = make_elements( xvec, nv );
00155                 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00156                 break;
00157             case 4:  // Read in tally data, make, and tag elements
00158                 mode++;
00159                 moab::EntityHandle elemhandle;
00160 
00161                 moab::EntityHandle vstart, vijk;
00162                 moab::EntityHandle connect[8];
00163                 // double d[3];
00164 
00165                 // vstart = MCNP_vertices.front();
00166                 vstart = *( vert_handles.begin() );
00167 
00168                 for( int i = 0; i < nv[0] - 1; i++ )
00169                 {
00170                     for( int j = 0; j < nv[1] - 1; j++ )
00171                     {
00172                         for( int k = 0; k < nv[2] - 1; k++ )
00173                         {
00174                             vijk = vstart + ( i + j * nv[0] + k * nv[0] * nv[1] );
00175 
00176                             // std::cout << vijk << std::endl;
00177 
00178                             connect[0] = vijk;
00179                             connect[1] = vijk + 1;
00180                             connect[2] = vijk + 1 + nv[0];
00181                             connect[3] = vijk + nv[0];
00182                             connect[4] = vijk + nv[0] * nv[1];
00183                             connect[5] = vijk + 1 + nv[0] * nv[1];
00184                             connect[6] = vijk + 1 + nv[0] + nv[0] * nv[1];
00185                             connect[7] = vijk + nv[0] + nv[0] * nv[1];
00186 
00187                             MBresult = MBI->create_element( moab::MBHEX, connect, 8, elemhandle );
00188                             if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
00189                             elem_handles.insert( elemhandle );
00190 
00191                             mcnpfile.getline( line, 10000 );
00192                             result = extract_tally_data( line, elemhandle );
00193                             if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00194                         }
00195                     }
00196                 }
00197 
00198                 /*
00199                                   for (MBRange::iterator rit=vert_handles.begin(); rit !=
00200                    vert_handles.end(); ++rit) { std::cout << *rit << std::endl;
00201                                   }
00202 
00203 
00204                                   for (int i=0; i < nv[0]-1; i++) {
00205                                     for (int j=0; j < nv[1]-1; j++) {
00206                                       for (int k=0; k < nv[2]-1; k++) {
00207                                         vijk = vstart + (i + j*nv[0] + k*nv[0]*nv[1]);
00208                                         connect[0] = vijk;
00209                                         connect[1] = vijk + 1;
00210                                         connect[2] = vijk + 1 + nv[0];
00211                                         connect[3] = vijk + nv[0];
00212                                         connect[4] = vijk + nv[0]*nv[1];
00213                                         connect[5] = vijk + 1 + nv[0]*nv[1];
00214                                         connect[6] = vijk + 1 + nv[0] + nv[0]*nv[1];
00215                                         connect[7] = vijk + nv[0] + nv[0]*nv[1];
00216 
00217                                         MBresult = MBI->create_element(MBHEX, connect, 8,
00218                    elemhandle); if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
00219                                         elem_handles.insert(elemhandle);
00220 
00221                                         mcnpfile.getline(line, 10000);
00222                                         result = extract_tally_data(line, elemhandle);
00223                                         if (result == MCNP_FAILURE) return MCNP_FAILURE;
00224 
00225                                     }
00226                                   }
00227                                 }
00228                 */
00229                 break;
00230             case 5:  // Ckeck for weirdness at end of file
00231                 if( !mcnpfile.eof() ) return MCNP_FAILURE;
00232                 break;
00233         }
00234     }
00235 
00236     std::cout << "SUCCESS! Read in " << elem_handles.size() << " elements!" << std::endl << std::endl;
00237     // MCNP_vertices.clear();
00238     vert_handles.clear();
00239     MCNP_elems.clear();
00240     return MCNP_SUCCESS;
00241 }
00242 
00243 MCNPError McnpData::read_coord_system( std::string s )
00244 {
00245 
00246     if( ( s.find( "Box" ) < 100 ) || ( s.find( "xyz" ) < 100 ) )
00247         coord_system = CARTESIAN;
00248     else if( s.find( "Cyl" ) < 100 )
00249         coord_system = CYLINDRICAL;
00250     else if( s.find( "Sph" ) < 100 )
00251         coord_system = SPHERICAL;
00252     else
00253         return MCNP_FAILURE;
00254 
00255     return MCNP_SUCCESS;
00256 }
00257 
00258 MCNPError McnpData::read_rotation_matrix( std::string s, int i )
00259 {
00260 
00261     int fpos = 0;
00262     MCNPError result;
00263 
00264     for( int j = 0; j < 4; j++ )
00265     {
00266         result = next_number( s, rotation_matrix[4 * i + j], fpos );
00267         if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00268     }
00269 
00270     return MCNP_SUCCESS;
00271 }
00272 
00273 MCNPError McnpData::make_elements( std::vector< double > x[3], int* n )
00274 {
00275 
00276     // double v[3];
00277     // MBEntityHandle dumhandle;
00278     // MBEntityHandle vstart, vijk;
00279     unsigned int num_verts = n[0] * n[1] * n[2];
00280     double* coords;
00281     coords = new double[3 * num_verts];
00282 
00283     /*
00284           // Enter the vertices ...
00285           for (int k=0; k < n[2]; k++) {
00286                 v[2] = x[2].at(k);
00287                 for (int j=0; j < n[1]; j++) {
00288                       v[1] = x[1].at(j);
00289                       for (int i=0; i < n[0]; i++) {
00290                             v[0] = x[0].at(i);
00291                             MBresult = MBI->create_vertex(v, dumhandle);
00292                             if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
00293                             MCNP_vertices.push_back(dumhandle);
00294 
00295                       }
00296                 }
00297           }
00298     */
00299 
00300     // Enter the vertices ...
00301     for( int k = 0; k < n[2]; k++ )
00302     {
00303         for( int j = 0; j < n[1]; j++ )
00304         {
00305             for( int i = 0; i < n[0]; i++ )
00306             {
00307                 unsigned int ijk = 3 * ( k * n[0] * n[1] + j * n[0] + i );
00308                 coords[ijk]      = x[0][i];
00309                 coords[ijk + 1]  = x[1][j];
00310                 coords[ijk + 2]  = x[2][k];
00311 
00312                 // std::cout << coords[ijk] << " " << coords[ijk+1] << " "
00313                 //          << coords[ijk+2] << std::endl;
00314 
00315                 // MCNP_vertices.push_back(dumhandle);
00316             }
00317         }
00318     }
00319 
00320     MBI->create_vertices( coords, num_verts, vert_handles );
00321 
00322     delete[] coords;
00323     return MCNP_SUCCESS;
00324 }
00325 
00326 MCNPError McnpData::initialize_tags()
00327 {
00328 
00329     MBI->tag_get_handle( TALLY_TAG, 1, moab::MB_TYPE_DOUBLE, tally_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
00330     MBI->tag_get_handle( ERROR_TAG, 1, moab::MB_TYPE_DOUBLE, relerr_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
00331 
00332     return MCNP_SUCCESS;
00333 }
00334 
00335 MCNPError McnpData::extract_tally_data( std::string s, moab::EntityHandle handle )
00336 {
00337 
00338     int fpos = 0;
00339     double d = 0;
00340 
00341     MCNPError result;
00342     moab::ErrorCode MBresult;
00343 
00344     // Discard first three lines
00345     for( int i = 0; i < 3; i++ )
00346     {
00347         result = next_number( s, d, fpos );
00348         if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00349     }
00350     // Need to read in tally entry and tag ...
00351     result = next_number( s, d, fpos );
00352     if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00353     MBresult = MBI->tag_set_data( tally_tag, &handle, 1, &d );
00354     if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
00355 
00356     // Need to read in relative error entry and tag ...
00357     result = next_number( s, d, fpos );
00358     if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00359     MBresult = MBI->tag_set_data( relerr_tag, &handle, 1, &d );
00360     if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
00361 
00362     return MCNP_SUCCESS;
00363 }
00364 
00365 MCNPError next_number( std::string s, double& d, int& p )
00366 {
00367 
00368     unsigned int slen = s.length();
00369     unsigned int j;
00370 
00371     for( unsigned int i = p; i < slen; i++ )
00372     {
00373         if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
00374         {
00375 
00376             j = s.find( " ", i );
00377 
00378             if( j > slen ) j = slen;
00379 
00380             // Extract the number out of the string
00381             d = std::atof( s.substr( i, j - i ).c_str() );
00382             p = j + 1;
00383 
00384             return MCNP_SUCCESS;
00385         }
00386     }
00387 
00388     return DONE;
00389 }
00390 
00391 int how_many_numbers( std::string s )
00392 {
00393 
00394     int n            = -1;
00395     int fpos         = 0;
00396     double d         = 0;
00397     MCNPError result = MCNP_SUCCESS;
00398 
00399     while( result != DONE )
00400     {
00401         result = next_number( s, d, fpos );
00402         n++;
00403     }
00404 
00405     return n;
00406 }
00407 
00408 MCNPError read_numbers( std::string s, int n, std::vector< double >& x )
00409 {
00410 
00411     MCNPError result;
00412     int fpos = 0;
00413     double d;
00414 
00415     for( int i = 0; i < n; i++ )
00416     {
00417         result = next_number( s, d, fpos );
00418         if( result == MCNP_FAILURE ) return MCNP_FAILURE;
00419         x.push_back( d );
00420     }
00421 
00422     return MCNP_SUCCESS;
00423 }
00424 
00425 MCNPError McnpData::transform_point( double* p, double* r, int csys, double* rmat )
00426 {
00427 
00428     double q[3];
00429 
00430     // Apply the rotation matrix
00431     for( unsigned int i = 0; i < 3; i++ )
00432     {
00433         q[i] = p[0] * rmat[4 * i] + p[1] * rmat[4 * i + 1] + p[2] * rmat[4 * i + 2] + rmat[4 * i + 3];
00434     }
00435 
00436     // Transform coordinate system
00437     switch( csys )
00438     {
00439         case CARTESIAN:
00440             r[0] = q[0];
00441             r[1] = q[1];
00442             r[2] = q[2];  // x, y, z
00443             break;
00444         case CYLINDRICAL:
00445             r[0] = sqrt( q[0] * q[0] + q[1] * q[1] );  // r
00446             r[1] = q[2];                               // z
00447             r[2] = c2pi * ( atan2( q[1], q[0] ) );     // theta (in rotations)
00448             break;
00449         case SPHERICAL:
00450             return MCNP_FAILURE;
00451         // break;
00452         default:
00453             return MCNP_FAILURE;
00454             // break;
00455     }
00456 
00457     return MCNP_SUCCESS;
00458 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines