mcnpmit.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <cstdlib>
#include "mcnpmit.hpp"
#include "moab/CartVect.hpp"
#include <cmath>
 
moab::Interface* mb_instance();
 
// Parameters
// const double pi = 3.141592653589793;
const double c2pi = 0.1591549430918954;
// const double cpi = 0.3183098861837907;
 
MCNPError next_number( std::string, double&, int& );
int how_many_numbers( std::string );
MCNPError read_numbers( std::string, int, std::vector< double >& );
 
// Constructor
McnpData::McnpData()
{
 
 // Default value for coordinate system
 coord_system = 0;
 
 // Default rotation matrix is identity matrix
 for( int i = 0; i < 4; i++ )
 {
 for( int j = 0; j < 4; j++ )
 {
 if( i == j )
 rotation_matrix[4 * i + j] = 1;
 else
 rotation_matrix[4 * i + j] = 0;
 }
 }
}
 
// Destructor
McnpData::~McnpData()
{
 
 // Vertices and elements
 MCNP_vertices.clear();
}
 
// Setting and retrieving coordinate sysem
MCNPError McnpData::set_coord_system( int k )
{
 coord_system = k;
 return MCNP_SUCCESS;
}
int McnpData::get_coord_system()
{
 return coord_system;
}
 
// Setting and retrieving roation matrix
MCNPError McnpData::set_rotation_matrix( double r[16] )
{
 for( int i = 0; i < 16; i++ )
 {
 rotation_matrix[i] = r[i];
 }
 return MCNP_SUCCESS;
}
double* McnpData::get_rotation_matrix()
{
 return rotation_matrix;
}
 
// Set the filename
MCNPError McnpData::set_filename( std::string fname )
{
 MCNP_filename = fname;
 return MCNP_SUCCESS;
}
std::string McnpData::get_filename()
{
 return MCNP_filename;
}
 
// Reading the MCNP file
MCNPError McnpData::read_mcnpfile( bool skip_mesh )
{
 
 MCNPError result;
 moab::ErrorCode MBresult;
 moab::CartVect tvect;
 
 std::vector< double > xvec[3];
 
 // Open the file
 std::ifstream mcnpfile;
 mcnpfile.open( MCNP_filename.c_str() );
 if( !mcnpfile )
 {
 std::cout << "Unable to open MCNP data file." << std::endl;
 return MCNP_FAILURE;
 }
 std::cout << std::endl;
 std::cout << "Reading MCNP input file..." << std::endl;
 
 // Prepare for file reading ...
 char line[10000];
 int mode = 0; // Set the file reading mode to read proper data
 int nv[3];
 
 // Read in the file ...
 while( !mcnpfile.eof() )
 {
 
 mcnpfile.getline( line, 10000 );
 // std::cout << line << std::endl;
 
 switch( mode )
 {
 case 0: // First line is a title
 mode++;
 break;
 case 1: // Coordinate system
 mode++;
 result = read_coord_system( line );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 break;
 case 2: // Rotation matrix
 mode++;
 for( int i = 0; i < 4; i++ )
 {
 mcnpfile.getline( line, 10000 );
 result = read_rotation_matrix( line, i );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 }
 if( skip_mesh ) return MCNP_SUCCESS;
 break;
 case 3: // Read in vertices and build elements
 mode++;
 
 for( int i = 0; i < 3; i++ )
 {
 // How many points in the x[i]-direction
 nv[i] = how_many_numbers( line );
 if( nv[i] <= 0 ) return MCNP_FAILURE;
 
 // Get space and read in these points
 result = read_numbers( line, nv[i], xvec[i] );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 
 // Update to the next line
 mcnpfile.getline( line, 10000 );
 }
 
 // Make the elements and vertices
 result = make_elements( xvec, nv );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 break;
 case 4: // Read in tally data, make, and tag elements
 mode++;
 moab::EntityHandle elemhandle;
 
 moab::EntityHandle vstart, vijk;
 moab::EntityHandle connect[8];
 // double d[3];
 
 // vstart = MCNP_vertices.front();
 vstart = *( vert_handles.begin() );
 
 for( int i = 0; i < nv[0] - 1; i++ )
 {
 for( int j = 0; j < nv[1] - 1; j++ )
 {
 for( int k = 0; k < nv[2] - 1; k++ )
 {
 vijk = vstart + ( i + j * nv[0] + k * nv[0] * nv[1] );
 
 // std::cout << vijk << std::endl;
 
 connect[0] = vijk;
 connect[1] = vijk + 1;
 connect[2] = vijk + 1 + nv[0];
 connect[3] = vijk + nv[0];
 connect[4] = vijk + nv[0] * nv[1];
 connect[5] = vijk + 1 + nv[0] * nv[1];
 connect[6] = vijk + 1 + nv[0] + nv[0] * nv[1];
 connect[7] = vijk + nv[0] + nv[0] * nv[1];
 
 MBresult = MBI->create_element( moab::MBHEX, connect, 8, elemhandle );
 if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
 elem_handles.insert( elemhandle );
 
 mcnpfile.getline( line, 10000 );
 result = extract_tally_data( line, elemhandle );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 }
 }
 }
 
 /*
 for (MBRange::iterator rit=vert_handles.begin(); rit !=
 vert_handles.end(); ++rit) { std::cout << *rit << std::endl;
 }
 
 
 for (int i=0; i < nv[0]-1; i++) {
 for (int j=0; j < nv[1]-1; j++) {
 for (int k=0; k < nv[2]-1; k++) {
 vijk = vstart + (i + j*nv[0] + k*nv[0]*nv[1]);
 connect[0] = vijk;
 connect[1] = vijk + 1;
 connect[2] = vijk + 1 + nv[0];
 connect[3] = vijk + nv[0];
 connect[4] = vijk + nv[0]*nv[1];
 connect[5] = vijk + 1 + nv[0]*nv[1];
 connect[6] = vijk + 1 + nv[0] + nv[0]*nv[1];
 connect[7] = vijk + nv[0] + nv[0]*nv[1];
 
 MBresult = MBI->create_element(MBHEX, connect, 8,
 elemhandle); if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
 elem_handles.insert(elemhandle);
 
 mcnpfile.getline(line, 10000);
 result = extract_tally_data(line, elemhandle);
 if (result == MCNP_FAILURE) return MCNP_FAILURE;
 
 }
 }
 }
 */
 break;
 case 5: // Ckeck for weirdness at end of file
 if( !mcnpfile.eof() ) return MCNP_FAILURE;
 break;
 }
 }
 
 std::cout << "SUCCESS! Read in " << elem_handles.size() << " elements!" << std::endl << std::endl;
 // MCNP_vertices.clear();
 vert_handles.clear();
 MCNP_elems.clear();
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::read_coord_system( std::string s )
{
 
 if( ( s.find( "Box" ) < 100 ) || ( s.find( "xyz" ) < 100 ) )
 coord_system = CARTESIAN;
 else if( s.find( "Cyl" ) < 100 )
 coord_system = CYLINDRICAL;
 else if( s.find( "Sph" ) < 100 )
 coord_system = SPHERICAL;
 else
 return MCNP_FAILURE;
 
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::read_rotation_matrix( std::string s, int i )
{
 
 int fpos = 0;
 MCNPError result;
 
 for( int j = 0; j < 4; j++ )
 {
 result = next_number( s, rotation_matrix[4 * i + j], fpos );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 }
 
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::make_elements( std::vector< double > x[3], int* n )
{
 
 // double v[3];
 // MBEntityHandle dumhandle;
 // MBEntityHandle vstart, vijk;
 unsigned int num_verts = n[0] * n[1] * n[2];
 double* coords;
 coords = new double[3 * num_verts];
 
 /*
 // Enter the vertices ...
 for (int k=0; k < n[2]; k++) {
 v[2] = x[2].at(k);
 for (int j=0; j < n[1]; j++) {
 v[1] = x[1].at(j);
 for (int i=0; i < n[0]; i++) {
 v[0] = x[0].at(i);
 MBresult = MBI->create_vertex(v, dumhandle);
 if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
 MCNP_vertices.push_back(dumhandle);
 
 }
 }
 }
 */
 
 // Enter the vertices ...
 for( int k = 0; k < n[2]; k++ )
 {
 for( int j = 0; j < n[1]; j++ )
 {
 for( int i = 0; i < n[0]; i++ )
 {
 unsigned int ijk = 3 * ( k * n[0] * n[1] + j * n[0] + i );
 coords[ijk] = x[0][i];
 coords[ijk + 1] = x[1][j];
 coords[ijk + 2] = x[2][k];
 
 // std::cout << coords[ijk] << " " << coords[ijk+1] << " "
 // << coords[ijk+2] << std::endl;
 
 // MCNP_vertices.push_back(dumhandle);
 }
 }
 }
 
 MBI->create_vertices( coords, num_verts, vert_handles );
 
 delete[] coords;
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::initialize_tags()
{
 
 MBI->tag_get_handle( TALLY_TAG, 1, moab::MB_TYPE_DOUBLE, tally_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
 MBI->tag_get_handle( ERROR_TAG, 1, moab::MB_TYPE_DOUBLE, relerr_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
 
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::extract_tally_data( std::string s, moab::EntityHandle handle )
{
 
 int fpos = 0;
 double d = 0;
 
 MCNPError result;
 moab::ErrorCode MBresult;
 
 // Discard first three lines
 for( int i = 0; i < 3; i++ )
 {
 result = next_number( s, d, fpos );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 }
 // Need to read in tally entry and tag ...
 result = next_number( s, d, fpos );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 MBresult = MBI->tag_set_data( tally_tag, &handle, 1, &d );
 if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
 
 // Need to read in relative error entry and tag ...
 result = next_number( s, d, fpos );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 MBresult = MBI->tag_set_data( relerr_tag, &handle, 1, &d );
 if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
 
 return MCNP_SUCCESS;
}
 
MCNPError next_number( std::string s, double& d, int& p )
{
 
 unsigned int slen = s.length();
 unsigned int j;
 
 for( unsigned int i = p; i < slen; i++ )
 {
 if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
 {
 
 j = s.find( " ", i );
 
 if( j > slen ) j = slen;
 
 // Extract the number out of the string
 d = std::atof( s.substr( i, j - i ).c_str() );
 p = j + 1;
 
 return MCNP_SUCCESS;
 }
 }
 
 return DONE;
}
 
int how_many_numbers( std::string s )
{
 
 int n = -1;
 int fpos = 0;
 double d = 0;
 MCNPError result = MCNP_SUCCESS;
 
 while( result != DONE )
 {
 result = next_number( s, d, fpos );
 n++;
 }
 
 return n;
}
 
MCNPError read_numbers( std::string s, int n, std::vector< double >& x )
{
 
 MCNPError result;
 int fpos = 0;
 double d;
 
 for( int i = 0; i < n; i++ )
 {
 result = next_number( s, d, fpos );
 if( result == MCNP_FAILURE ) return MCNP_FAILURE;
 x.push_back( d );
 }
 
 return MCNP_SUCCESS;
}
 
MCNPError McnpData::transform_point( double* p, double* r, int csys, double* rmat )
{
 
 double q[3];
 
 // Apply the rotation matrix
 for( unsigned int i = 0; i < 3; i++ )
 {
 q[i] = p[0] * rmat[4 * i] + p[1] * rmat[4 * i + 1] + p[2] * rmat[4 * i + 2] + rmat[4 * i + 3];
 }
 
 // Transform coordinate system
 switch( csys )
 {
 case CARTESIAN:
 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
 r[1] = q[2]; // z
 r[2] = c2pi * ( atan2( q[1], q[0] ) ); // theta (in rotations)
 break;
 case SPHERICAL:
 return MCNP_FAILURE;
 // break;
 default:
 return MCNP_FAILURE;
 // break;
 }
 
 return MCNP_SUCCESS;
}