moab::ReadMCNP5 Class Reference

#include <ReadMCNP5.hpp>

Inheritance diagram for moab::ReadMCNP5:

Collaboration diagram for moab::ReadMCNP5:

Public Member Functions
ErrorCode	load_file (const char *file_name, const EntityHandle *file_set, const FileOptions &opts, const SubsetList *subset_list=0, const Tag *file_id_tag=0)
	Load mesh from a file. More...
ErrorCode	read_tag_values (const char *file_name, const char *tag_name, const FileOptions &opts, std::vector< int > &tag_values_out, const SubsetList *subset_list=0)
	Read tag values from a file. More...
	ReadMCNP5 (Interface *impl=NULL)
virtual	~ReadMCNP5 ()
Public Member Functions inherited from moab::ReaderIface
virtual	~ReaderIface ()

Static Public Member Functions
static ReaderIface *	factory (Interface *)

Private Types
enum	coordinate_system { NO_SYSTEM , CARTESIAN , CYLINDRICAL , SPHERICAL }
enum	particle { NEUTRON , PHOTON , ELECTRON }

Private Member Functions
ErrorCode	load_one_file (const char *fname, const EntityHandle *input_meshset, const FileOptions &options, const bool average)
ErrorCode	create_tags (Tag &date_and_time_tag, Tag &title_tag, Tag &nps_tag, Tag &tally_number_tag, Tag &tally_comment_tag, Tag &tally_particle_tag, Tag &tally_coord_sys_tag, Tag &tally_tag, Tag &error_tag)
ErrorCode	read_file_header (std::fstream &file, bool debug, char date_and_time[100], char title[100], unsigned long int &nps)
ErrorCode	set_header_tags (EntityHandle output_meshset, char date_and_time[100], char title[100], unsigned long int nps, Tag data_and_time_tag, Tag title_tag, Tag nps_tag)
ErrorCode	read_tally_header (std::fstream &file, bool debug, unsigned int &tally_number, char tally_comment[100], particle &tally_particle)
ErrorCode	get_tally_particle (std::string a, bool debug, particle &tally_particle)
ErrorCode	read_mesh_planes (std::fstream &file, bool debug, std::vector< double > planes[3], coordinate_system &coord_sys)
ErrorCode	get_mesh_plane (std::istringstream &ss, bool debug, std::vector< double > &plane)
ErrorCode	read_element_values_and_errors (std::fstream &file, bool debug, std::vector< double > planes[3], unsigned int n_chopped_x0_planes, unsigned int n_chopped_x2_planes, particle tally_particle, double values[], double errors[])
ErrorCode	set_tally_tags (EntityHandle tally_meshset, unsigned int tally_number, char tally_comment[100], particle tally_particle, coordinate_system tally_coord_sys, Tag tally_number_tag, Tag tally_comment_tag, Tag tally_particle_tag, Tag tally_coord_sys_tag)
ErrorCode	create_vertices (std::vector< double > planes[3], bool debug, EntityHandle &start_vert, coordinate_system coord_sys, EntityHandle tally_meshset)
ErrorCode	create_elements (bool debug, std::vector< double > planes[3], unsigned int n_chopped_x0_planes, unsigned int n_chopped_x2_planes, EntityHandle start_vert, double values[], double errors[], Tag tally_tag, Tag error_tag, EntityHandle tally_meshset, coordinate_system tally_coord_sys)
ErrorCode	average_with_existing_tally (bool debug, unsigned long int &new_nps, unsigned long int nps, unsigned int tally_number, Tag tally_number_tag, Tag nps_tag, Tag tally_tag, Tag error_tag, double values[], double errors[], unsigned int n_elements)
ErrorCode	transform_point_to_cartesian (double *in, double *out, coordinate_system coord_sys)
ErrorCode	average_tally_values (const unsigned long int nps0, const unsigned long int nps1, double *values0, const double *values1, double *errors0, const double *errors1, const unsigned long int n_values)

Private Attributes
ReadUtilIface *	readMeshIface
Interface *	MBI
const Tag *	fileIDTag
int	nodeId
int	elemId

Static Private Attributes
static const double	PI = 3.141592653589793
static const double	C2PI = 0.1591549430918954
static const double	CPI = 0.3183098861837907

Detailed Description

Definition at line 52 of file ReadMCNP5.hpp.

Member Enumeration Documentation

coordinate_system

enum moab::ReadMCNP5::coordinate_system

private

Enumerator
NO_SYSTEM
CARTESIAN
CYLINDRICAL
SPHERICAL

Definition at line 84 of file ReadMCNP5.hpp.

 {
 NO_SYSTEM,
 CARTESIAN,
 CYLINDRICAL,
 SPHERICAL
 };

particle

enum moab::ReadMCNP5::particle

private

Enumerator
NEUTRON
PHOTON
ELECTRON

Definition at line 91 of file ReadMCNP5.hpp.

 {
 NEUTRON,
 PHOTON,
 ELECTRON
 };

Constructor & Destructor Documentation

ReadMCNP5()

moab::ReadMCNP5::ReadMCNP5

(

Interface *

impl = NULL

)

Definition at line 46 of file ReadMCNP5.cpp.

 : MBI( impl ), fileIDTag( NULL ), nodeId( 0 ), elemId( 0 )
{
 assert( NULL != impl );
 MBI->query_interface( readMeshIface );
 assert( NULL != readMeshIface );
}

References MBI, moab::Interface::query_interface(), and readMeshIface.

Referenced by factory().

~ReadMCNP5()

moab::ReadMCNP5::~ReadMCNP5 ( )

virtual

Definition at line 54 of file ReadMCNP5.cpp.

{
 if( readMeshIface )
 {
 MBI->release_interface( readMeshIface );
 readMeshIface = 0;
 }
}

References MBI, readMeshIface, and moab::Interface::release_interface().

Member Function Documentation

average_tally_values()

ErrorCode moab::ReadMCNP5::average_tally_values ( const unsigned long int  nps0, const unsigned long int  nps1, double *  values0, const double *  values1, double *  errors0, const double *  errors1, const unsigned long int  n_values  )

private

Definition at line 1006 of file ReadMCNP5.cpp.

{
 for( unsigned long int i = 0; i < n_values; i++ )
 {
 // std::cout << " values0=" << values0[i] << " values1=" << values1[i]
 // << " errors0=" << errors0[i] << " errors1=" << errors1[i]
 // << " nps0=" << nps0 << " nps1=" << nps1 << std::endl;
 errors0[i] = sqrt( pow( values0[i] * errors0[i] * nps0, 2 ) + pow( values1[i] * errors1[i] * nps1, 2 ) ) /
 ( values0[i] * nps0 + values1[i] * nps1 );
 
 // It is possible to introduce nans if the values are zero.
 if( !Util::is_finite( errors0[i] ) ) errors0[i] = 1.0;
 
 values0[i] = ( values0[i] * nps0 + values1[i] * nps1 ) / ( nps0 + nps1 );
 
 // std::cout << " values0=" << values0[i] << " errors0=" << errors0[i] << std::endl;
 }
 // REMEMBER TO UPDATE NPS0 = NPS0 + NPS1 after this
 
 return MB_SUCCESS;
}

References moab::Util::is_finite(), and MB_SUCCESS.

Referenced by average_with_existing_tally().

average_with_existing_tally()

ErrorCode moab::ReadMCNP5::average_with_existing_tally ( bool  debug, unsigned long int &  new_nps, unsigned long int  nps, unsigned int  tally_number, Tag  tally_number_tag, Tag  nps_tag, Tag  tally_tag, Tag  error_tag, double  values[], double  errors[], unsigned int  n_elements  )

private

Definition at line 883 of file ReadMCNP5.cpp.

{
 // Get the tally number
 ErrorCode result;
 
 // Match the tally number with one from the existing meshtal file
 Range matching_tally_number_sets;
 const void* const tally_number_val[] = { &tally_number };
 result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &tally_number_tag, tally_number_val, 1,
 matching_tally_number_sets );
 if( MB_SUCCESS != result ) return result;
 if( debug ) std::cout << "number of matching meshsets=" << matching_tally_number_sets.size() << std::endl;
 assert( 1 == matching_tally_number_sets.size() );
 
 // Identify which of the meshsets is existing
 EntityHandle existing_meshset;
 existing_meshset = matching_tally_number_sets.front();
 
 // Get the existing elements from the set
 Range existing_elements;
 result = MBI->get_entities_by_type( existing_meshset, MBHEX, existing_elements );
 if( MB_SUCCESS != result ) return result;
 assert( existing_elements.size() == n_elements );
 
 // Get the nps of the existing and new tally
 unsigned long int nps0;
 Range sets_with_this_tag;
 result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, sets_with_this_tag );
 if( MB_SUCCESS != result ) return result;
 if( debug ) std::cout << "number of nps sets=" << sets_with_this_tag.size() << std::endl;
 assert( 1 == sets_with_this_tag.size() );
 result = MBI->tag_get_data( nps_tag, &sets_with_this_tag.front(), 1, &nps0 );
 if( MB_SUCCESS != result ) return result;
 if( debug ) std::cout << "nps0=" << nps0 << " nps1=" << nps1 << std::endl;
 new_nps = nps0 + nps1;
 
 // Get tally values from the existing elements
 double* values0 = new double[existing_elements.size()];
 double* errors0 = new double[existing_elements.size()];
 result = MBI->tag_get_data( tally_tag, existing_elements, values0 );
 if( MB_SUCCESS != result )
 {
 delete[] values0;
 delete[] errors0;
 return result;
 }
 result = MBI->tag_get_data( error_tag, existing_elements, errors0 );
 if( MB_SUCCESS != result )
 {
 delete[] values0;
 delete[] errors0;
 return result;
 }
 
 // Average the values and errors
 result = average_tally_values( nps0, nps1, values0, values1, errors0, errors1, n_elements );
 if( MB_SUCCESS != result )
 {
 delete[] values0;
 delete[] errors0;
 return result;
 }
 
 // Set the averaged information back onto the existing elements
 result = MBI->tag_set_data( tally_tag, existing_elements, values0 );
 if( MB_SUCCESS != result )
 {
 delete[] values0;
 delete[] errors0;
 return result;
 }
 result = MBI->tag_set_data( error_tag, existing_elements, errors0 );
 if( MB_SUCCESS != result )
 {
 delete[] values0;
 delete[] errors0;
 return result;
 }
 
 // Cleanup
 delete[] values0;
 delete[] errors0;
 
 return MB_SUCCESS;
}

References average_tally_values(), moab::debug, ErrorCode, moab::Range::front(), moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), MB_SUCCESS, MBENTITYSET, MBHEX, MBI, moab::Range::size(), moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().

Referenced by load_one_file().

create_elements()

ErrorCode moab::ReadMCNP5::create_elements ( bool  debug, std::vector< double >  planes[3], unsigned int  n_chopped_x0_planes, unsigned int  n_chopped_x2_planes, EntityHandle  start_vert, double  values[], double  errors[], Tag  tally_tag, Tag  error_tag, EntityHandle  tally_meshset, coordinate_system  tally_coord_sys  )

private

Definition at line 793 of file ReadMCNP5.cpp.

{
 ErrorCode result;
 unsigned int index;
 EntityHandle start_element = 0;
 unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
 EntityHandle* connect;
 result = readMeshIface->get_element_connect( n_elements, 8, MBHEX, MB_START_ID, start_element, connect );
 if( MB_SUCCESS != result ) return result;
 assert( 0 != start_element ); // Check for NULL
 
 unsigned int counter = 0;
 for( unsigned int i = 0; i < planes[0].size() - 1; i++ )
 {
 for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
 {
 for( unsigned int k = 0; k < planes[2].size() - 1; k++ )
 {
 index = start_vert + i + j * planes[0].size() + k * planes[0].size() * planes[1].size();
 // For rectangular mesh, the file prints: x y z
 // z changes the fastest and x changes the slowest.
 // This means that the connectivity ordering is not consistent between
 // rectangular and cylindrical mesh.
 if( CARTESIAN == tally_coord_sys )
 {
 connect[0] = index;
 connect[1] = index + 1;
 connect[2] = index + 1 + planes[0].size();
 connect[3] = index + planes[0].size();
 connect[4] = index + planes[0].size() * planes[1].size();
 connect[5] = index + 1 + planes[0].size() * planes[1].size();
 connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
 connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
 // For cylindrical mesh, the file prints: r z theta
 // Theta changes the fastest and r changes the slowest.
 }
 else if( CYLINDRICAL == tally_coord_sys )
 {
 connect[0] = index;
 connect[1] = index + 1;
 connect[2] = index + 1 + planes[0].size() * planes[1].size();
 connect[3] = index + planes[0].size() * planes[1].size();
 connect[4] = index + planes[0].size();
 connect[5] = index + 1 + planes[0].size();
 connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
 connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
 }
 else
 return MB_NOT_IMPLEMENTED;
 
 connect += 8;
 counter++;
 }
 }
 }
 if( counter != n_elements ) std::cout << "counter=" << counter << " n_elements=" << n_elements << std::endl;
 
 Range element_range( start_element, start_element + n_elements - 1 );
 result = MBI->tag_set_data( tally_tag, element_range, values );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( error_tag, element_range, errors );
 if( MB_SUCCESS != result ) return result;
 
 // Add the elements to the tally set
 result = MBI->add_entities( tally_meshset, element_range );
 if( MB_SUCCESS != result ) return result;
 if( debug ) std::cout << "Read " << n_elements << " elements from tally." << std::endl;
 
 if( fileIDTag )
 {
 result = readMeshIface->assign_ids( *fileIDTag, element_range, elemId );
 if( MB_SUCCESS != result ) return result;
 elemId += element_range.size();
 }
 
 return MB_SUCCESS;
}

References moab::Interface::add_entities(), moab::ReadUtilIface::assign_ids(), CARTESIAN, CYLINDRICAL, moab::debug, elemId, ErrorCode, fileIDTag, moab::ReadUtilIface::get_element_connect(), MB_NOT_IMPLEMENTED, MB_START_ID, MB_SUCCESS, MBHEX, MBI, readMeshIface, size, moab::Range::size(), and moab::Interface::tag_set_data().

Referenced by load_one_file().

create_tags()

ErrorCode moab::ReadMCNP5::create_tags ( Tag &  date_and_time_tag, Tag &  title_tag, Tag &  nps_tag, Tag &  tally_number_tag, Tag &  tally_comment_tag, Tag &  tally_particle_tag, Tag &  tally_coord_sys_tag, Tag &  tally_tag, Tag &  error_tag  )

private

Definition at line 335 of file ReadMCNP5.cpp.

{
 ErrorCode result;
 result = MBI->tag_get_handle( "DATE_AND_TIME_TAG", 100, MB_TYPE_OPAQUE, date_and_time_tag,
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "TITLE_TAG", 100, MB_TYPE_OPAQUE, title_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "NPS_TAG", sizeof( unsigned long int ), MB_TYPE_OPAQUE, nps_tag,
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result =
 MBI->tag_get_handle( "TALLY_NUMBER_TAG", 1, MB_TYPE_INTEGER, tally_number_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "TALLY_COMMENT_TAG", 100, MB_TYPE_OPAQUE, tally_comment_tag,
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "TALLY_PARTICLE_TAG", sizeof( particle ), MB_TYPE_OPAQUE, tally_particle_tag,
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "TALLY_COORD_SYS_TAG", sizeof( coordinate_system ), MB_TYPE_OPAQUE,
 tally_coord_sys_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "TALLY_TAG", 1, MB_TYPE_DOUBLE, tally_tag, MB_TAG_DENSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_get_handle( "ERROR_TAG", 1, MB_TYPE_DOUBLE, error_tag, MB_TAG_DENSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result ) return result;
 
 return MB_SUCCESS;
}

References ErrorCode, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_SPARSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, MBI, and moab::Interface::tag_get_handle().

Referenced by load_one_file().

create_vertices()

ErrorCode moab::ReadMCNP5::create_vertices ( std::vector< double >  planes[3], bool  debug, EntityHandle &  start_vert, coordinate_system  coord_sys, EntityHandle  tally_meshset  )

private

Definition at line 741 of file ReadMCNP5.cpp.

{
 // The only info needed to build elements is the mesh plane boundaries.
 ErrorCode result;
 int n_verts = planes[0].size() * planes[1].size() * planes[2].size();
 if( debug ) std::cout << "n_verts=" << n_verts << std::endl;
 std::vector< double* > coord_arrays( 3 );
 result = readMeshIface->get_node_coords( 3, n_verts, MB_START_ID, start_vert, coord_arrays );
 if( MB_SUCCESS != result ) return result;
 assert( 0 != start_vert ); // Check for NULL
 
 for( unsigned int k = 0; k < planes[2].size(); k++ )
 {
 for( unsigned int j = 0; j < planes[1].size(); j++ )
 {
 for( unsigned int i = 0; i < planes[0].size(); i++ )
 {
 unsigned int idx = ( k * planes[0].size() * planes[1].size() + j * planes[0].size() + i );
 double in[3], out[3];
 
 in[0] = planes[0][i];
 in[1] = planes[1][j];
 in[2] = planes[2][k];
 result = transform_point_to_cartesian( in, out, coord_sys );
 if( MB_SUCCESS != result ) return result;
 
 // Cppcheck warning (false positive): variable coord_arrays is assigned a value that
 // is never used
 coord_arrays[0][idx] = out[0];
 coord_arrays[1][idx] = out[1];
 coord_arrays[2][idx] = out[2];
 }
 }
 }
 Range vert_range( start_vert, start_vert + n_verts - 1 );
 result = MBI->add_entities( tally_meshset, vert_range );
 if( MB_SUCCESS != result ) return result;
 
 if( fileIDTag )
 {
 result = readMeshIface->assign_ids( *fileIDTag, vert_range, nodeId );
 if( MB_SUCCESS != result ) return result;
 nodeId += vert_range.size();
 }
 
 return MB_SUCCESS;
}

References moab::Interface::add_entities(), moab::ReadUtilIface::assign_ids(), moab::debug, ErrorCode, fileIDTag, moab::ReadUtilIface::get_node_coords(), MB_START_ID, MB_SUCCESS, MBI, nodeId, readMeshIface, moab::Range::size(), and transform_point_to_cartesian().

Referenced by load_one_file().

factory()

ReaderIface * moab::ReadMCNP5::factory ( Interface *  iface )

static

Definition at line 40 of file ReadMCNP5.cpp.

{
 return new ReadMCNP5( iface );
}

References iface, and ReadMCNP5().

Referenced by moab::ReaderWriterSet::ReaderWriterSet().

get_mesh_plane()

ErrorCode moab::ReadMCNP5::get_mesh_plane ( std::istringstream &  ss, bool  debug, std::vector< double > &  plane  )

private

Definition at line 647 of file ReadMCNP5.cpp.

{
 double value;
 plane.clear();
 while( !ss.eof() )
 {
 ss >> value;
 plane.push_back( value );
 if( debug ) std::cout << value << " ";
 }
 if( debug ) std::cout << std::endl;
 
 return MB_SUCCESS;
}

References moab::debug, and MB_SUCCESS.

Referenced by read_mesh_planes().

get_tally_particle()

ErrorCode moab::ReadMCNP5::get_tally_particle ( std::string  a, bool  debug, particle &  tally_particle  )

private

Definition at line 479 of file ReadMCNP5.cpp.

{
 if( std::string::npos != a.find( "This is a neutron mesh tally." ) )
 {
 tally_particle = NEUTRON;
 }
 else if( std::string::npos != a.find( "This is a photon mesh tally." ) )
 {
 tally_particle = PHOTON;
 }
 else if( std::string::npos != a.find( "This is an electron mesh tally." ) )
 {
 tally_particle = ELECTRON;
 }
 else
 return MB_FAILURE;
 
 if( debug ) std::cout << "tally_particle=| " << tally_particle << std::endl;
 
 return MB_SUCCESS;
}

References moab::debug, ELECTRON, MB_SUCCESS, NEUTRON, and PHOTON.

Referenced by read_tally_header().

load_file()

ErrorCode moab::ReadMCNP5::load_file ( const char *  file_name, const EntityHandle *  file_set, const FileOptions &  opts, const SubsetList *  subset_list = 0, const Tag *  file_id_tag = 0  )

virtual

Load mesh from a file.

Method all readers must provide to import a mesh.

Parameters

file_name	The file to read.
file_set	Optional pointer to entity set representing file. If this is not NULL, reader may optionally tag the pointed-to set with format-specific meta-data.
subset_list	An optional struct pointer specifying the tags identifying entity sets to be read.
file_id_tag	If specified, reader should store for each entity it reads, a unique integer ID for this tag.

Author

Jason Kraftcheck

Implements moab::ReaderIface.

Definition at line 73 of file ReadMCNP5.cpp.

{
 // At this time there is no support for reading a subset of the file
 if( subset_list )
 {
 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for meshtal" );
 }
 
 nodeId = elemId = 0;
 fileIDTag = file_id_tag;
 
 // Average several meshtal files if the AVERAGE_TALLY option is given.
 // In this case, the integer value is the number of files to average.
 // If averaging, the filename passed to load_file is assumed to be the
 // root filename. The files are indexed as "root_filename""index".meshtal.
 // Indices start with 1.
 int n_files;
 bool average = false;
 ErrorCode result;
 if( MB_SUCCESS == options.get_int_option( "AVERAGE_TALLY", n_files ) )
 {
 // Read the first file (but do not average -> can't average a single file)
 result = load_one_file( filename, input_meshset, options, average );
 if( MB_SUCCESS != result ) return result;
 
 // Get the root filename
 std::string root_filename( filename );
 int length = root_filename.length();
 root_filename.erase( length - sizeof( ".meshtal" ) );
 
 // Average the first file with the rest of the files
 average = true;
 for( int i = 2; i <= n_files; i++ )
 {
 std::stringstream index;
 index << i;
 std::string subsequent_filename = root_filename + index.str() + ".meshtal";
 result = load_one_file( subsequent_filename.c_str(), input_meshset, options, average );
 if( MB_SUCCESS != result ) return result;
 }
 
 // If not averaging, read a single file
 }
 else
 {
 result = load_one_file( filename, input_meshset, options, average );
 if( MB_SUCCESS != result ) return result;
 }
 
 return MB_SUCCESS;
}

References elemId, ErrorCode, fileIDTag, moab::FileOptions::get_int_option(), length(), load_one_file(), MB_SET_ERR, MB_SUCCESS, MB_UNSUPPORTED_OPERATION, and nodeId.

load_one_file()

ErrorCode moab::ReadMCNP5::load_one_file ( const char *  fname, const EntityHandle *  input_meshset, const FileOptions &  options, const bool  average  )

private

Definition at line 131 of file ReadMCNP5.cpp.

{
 bool debug = false;
 if( debug ) std::cout << "begin ReadMCNP5::load_one_file" << std::endl;
 
 ErrorCode result;
 std::fstream file;
 file.open( fname, std::fstream::in );
 char line[10000];
 
 // Create tags
 Tag date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag, tally_particle_tag,
 tally_coord_sys_tag, tally_tag, error_tag;
 
 result = create_tags( date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag,
 tally_particle_tag, tally_coord_sys_tag, tally_tag, error_tag );
 if( MB_SUCCESS != result ) return result;
 
 // ******************************************************************
 // This info exists only at the top of each meshtal file
 // ******************************************************************
 
 // Define characteristics of the entire file
 char date_and_time[100] = "";
 char title[100] = "";
 // This file's number of particles
 unsigned long int nps;
 // Sum of this file's and existing file's nps for averaging
 unsigned long int new_nps;
 
 // Read the file header
 result = read_file_header( file, debug, date_and_time, title, nps );
 if( MB_SUCCESS != result ) return result;
 
 // Blank line
 file.getline( line, 10000 );
 
 // Everything stored in the file being read will be in the input_meshset.
 // If this is being saved in MOAB, set header tags
 if( !average && 0 != input_meshset )
 {
 result = set_header_tags( *input_meshset, date_and_time, title, nps, date_and_time_tag, title_tag, nps_tag );
 if( MB_SUCCESS != result ) return result;
 }
 
 // ******************************************************************
 // This info is repeated for each tally in the meshtal file.
 // ******************************************************************
 
 // If averaging, nps will hold the sum of particles simulated in both tallies.
 while( !file.eof() )
 {
 // Define characteristics of this tally
 unsigned int tally_number;
 char tally_comment[100] = "";
 particle tally_particle;
 coordinate_system tally_coord_sys;
 std::vector< double > planes[3];
 unsigned int n_chopped_x0_planes;
 unsigned int n_chopped_x2_planes;
 
 // Read tally header
 result = read_tally_header( file, debug, tally_number, tally_comment, tally_particle );
 if( MB_SUCCESS != result ) return result;
 
 // Blank line
 file.getline( line, 10000 );
 std::string l = line;
 // if this string is present then skip the following blank line
 if( std::string::npos != l.find( "This mesh tally is modified by a dose response function." ) )
 {
 file.getline( line, 10000 );
 }
 
 // Read mesh planes
 result = read_mesh_planes( file, debug, planes, tally_coord_sys );
 if( MB_SUCCESS != result ) return result;
 
 // Get energy boundaries
 file.getline( line, 10000 );
 std::string a = line;
 if( debug ) std::cout << "Energy bin boundaries:=| " << a << std::endl;
 
 // Blank
 file.getline( line, 10000 );
 
 // Column headers
 file.getline( line, 10000 );
 
 // If using cylindrical mesh, it may be necessary to chop off the last theta element.
 // We chop off the last theta plane because the elements will be wrong and skew up
 // the tree building code. This is
 // because the hex elements are a linear approximation to the cylindrical elements.
 // Chopping off the last plane is problem-dependent, and due to MCNP5's mandate
 // that the cylindrical mesh must span 360 degrees.
 if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_LAST_AZIMUTHAL_PLANE" ) )
 {
 planes[2].pop_back();
 n_chopped_x2_planes = 1;
 if( debug ) std::cout << "remove last cylindrical plane option found" << std::endl;
 }
 else
 {
 n_chopped_x2_planes = 0;
 }
 
 // If using cylindrical mesh, it may be necessary to chop off the first radial element.
 // These elements extend from the axis and often do not cover areas of interest. For
 // example, if the mesh was meant to cover r=390-400, the first layer will go from
 // 0-390 and serve as incorrect source elements for interpolation.
 if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_FIRST_RADIAL_PLANE" ) )
 {
 std::vector< double >::iterator front = planes[0].begin();
 planes[0].erase( front );
 n_chopped_x0_planes = 1;
 if( debug ) std::cout << "remove first radial plane option found" << std::endl;
 }
 else
 {
 n_chopped_x0_planes = 0;
 }
 
 // Read the values and errors of each element from the file.
 // Do not read values that are chopped off.
 unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
 double* values = new double[n_elements];
 double* errors = new double[n_elements];
 result = read_element_values_and_errors( file, debug, planes, n_chopped_x0_planes, n_chopped_x2_planes,
 tally_particle, values, errors );
 if( MB_SUCCESS != result ) return result;
 
 // Blank line
 file.getline( line, 10000 );
 
 // ****************************************************************
 // This tally has been read. If it is not being averaged, build tags,
 // vertices and elements. If it is being averaged, average the data
 // with a tally already existing in the MOAB instance.
 // ****************************************************************
 if( !average )
 {
 EntityHandle tally_meshset;
 result = MBI->create_meshset( MESHSET_SET, tally_meshset );
 if( MB_SUCCESS != result ) return result;
 
 // Set tags on the tally
 result = set_tally_tags( tally_meshset, tally_number, tally_comment, tally_particle, tally_coord_sys,
 tally_number_tag, tally_comment_tag, tally_particle_tag, tally_coord_sys_tag );
 if( MB_SUCCESS != result ) return result;
 
 // The only info needed to build elements is the mesh plane boundaries.
 // Build vertices...
 EntityHandle start_vert = 0;
 result = create_vertices( planes, debug, start_vert, tally_coord_sys, tally_meshset );
 if( MB_SUCCESS != result ) return result;
 
 // Build elements and tag them with tally values and errors, then add
 // them to the tally_meshset.
 result = create_elements( debug, planes, n_chopped_x0_planes, n_chopped_x2_planes, start_vert, values,
 errors, tally_tag, error_tag, tally_meshset, tally_coord_sys );
 if( MB_SUCCESS != result ) return result;
 
 // Add this tally's meshset to the output meshset
 if( debug ) std::cout << "not averaging tally" << std::endl;
 
 // Average the tally values, then delete stuff that was created
 }
 else
 {
 if( debug ) std::cout << "averaging tally" << std::endl;
 result = average_with_existing_tally( debug, new_nps, nps, tally_number, tally_number_tag, nps_tag,
 tally_tag, error_tag, values, errors, n_elements );
 if( MB_SUCCESS != result ) return result;
 }
 
 // Clean up
 delete[] values;
 delete[] errors;
 }
 
 // If we are averaging, delete the remainder of this file's information.
 // Add the new nps to the existing file's nps if we are averaging.
 // This is calculated during every tally averaging but only used after the last one.
 if( average )
 {
 Range matching_nps_sets;
 result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, matching_nps_sets );
 if( MB_SUCCESS != result ) return result;
 if( debug ) std::cout << "number of matching nps meshsets=" << matching_nps_sets.size() << std::endl;
 assert( 1 == matching_nps_sets.size() );
 result = MBI->tag_set_data( nps_tag, matching_nps_sets, &new_nps );
 if( MB_SUCCESS != result ) return result;
 
 // If this file is not being averaged, return the output meshset.
 }
 
 file.close();
 return MB_SUCCESS;
}

References average_with_existing_tally(), create_elements(), moab::Interface::create_meshset(), create_tags(), create_vertices(), CYLINDRICAL, moab::debug, ErrorCode, moab::Interface::get_entities_by_type_and_tag(), moab::FileOptions::get_null_option(), MB_SUCCESS, MBENTITYSET, MBI, MESHSET_SET, read_element_values_and_errors(), read_file_header(), read_mesh_planes(), read_tally_header(), set_header_tags(), set_tally_tags(), size, moab::Range::size(), and moab::Interface::tag_set_data().

Referenced by load_file().

read_element_values_and_errors()

ErrorCode moab::ReadMCNP5::read_element_values_and_errors ( std::fstream &  file, bool  debug, std::vector< double >  planes[3], unsigned int  n_chopped_x0_planes, unsigned int  n_chopped_x2_planes, particle  tally_particle, double  values[], double  errors[]  )

private

Definition at line 662 of file ReadMCNP5.cpp.

{
 unsigned int index = 0;
 // Need to read every line in the file, even if we chop off some elements
 for( unsigned int i = 0; i < planes[0].size() - 1 + n_chopped_x0_planes; i++ )
 {
 for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
 {
 for( unsigned int k = 0; k < planes[2].size() - 1 + n_chopped_x2_planes; k++ )
 {
 char line[100];
 file.getline( line, 100 );
 // If this element has been chopped off, skip it
 if( i < n_chopped_x0_planes ) continue;
 if( k >= planes[2].size() - 1 && k < planes[2].size() - 1 + n_chopped_x2_planes ) continue;
 std::string a = line;
 std::stringstream ss( a );
 double centroid[3];
 double energy;
 // For some reason, photon tallies print the energy in the first column
 if( PHOTON == tally_particle ) ss >> energy;
 // The centroid is not used in this reader
 ss >> centroid[0];
 ss >> centroid[1];
 ss >> centroid[2];
 // Only the tally values and errors are used
 ss >> values[index];
 ss >> errors[index];
 
 // Make sure that input data is good
 if( !Util::is_finite( errors[index] ) )
 {
 std::cerr << "found nan error while reading file" << std::endl;
 errors[index] = 1.0;
 }
 if( !Util::is_finite( values[index] ) )
 {
 std::cerr << "found nan value while reading file" << std::endl;
 values[index] = 0.0;
 }
 
 index++;
 }
 }
 }
 
 return MB_SUCCESS;
}

References moab::Util::is_finite(), MB_SUCCESS, PHOTON, and size.

Referenced by load_one_file().

read_file_header()

ErrorCode moab::ReadMCNP5::read_file_header ( std::fstream &  file, bool  debug, char  date_and_time[100], char  title[100], unsigned long int &  nps  )

private

Definition at line 374 of file ReadMCNP5.cpp.

{
 // Get simulation date and time
 // mcnp version 5 ld=11242008 probid = 03/23/09 13:38:56
 char line[100];
 file.getline( line, 100 );
 date_and_time = line;
 if( debug ) std::cout << "date_and_time=| " << date_and_time << std::endl;
 
 // Get simulation title
 // iter Module 4
 file.getline( line, 100 );
 title = line;
 if( debug ) std::cout << "title=| " << title << std::endl;
 
 // Get number of histories
 // Number of histories used for normalizing tallies = 50000000.00
 file.getline( line, 100 );
 std::string a = line;
 std::string::size_type b = a.find( "Number of histories used for normalizing tallies =" );
 if( std::string::npos != b )
 {
 std::istringstream nps_ss(
 a.substr( b + sizeof( "Number of histories used for normalizing tallies =" ), 100 ) );
 nps_ss >> nps;
 if( debug ) std::cout << "nps=| " << nps << std::endl;
 }
 else
 return MB_FAILURE;
 
 return MB_SUCCESS;
}

References moab::debug, and MB_SUCCESS.

Referenced by load_one_file().

read_mesh_planes()

ErrorCode moab::ReadMCNP5::read_mesh_planes ( std::fstream &  file, bool  debug, std::vector< double >  planes[3], coordinate_system &  coord_sys  )

private

Definition at line 501 of file ReadMCNP5.cpp.

{
 // Tally bin boundaries:
 ErrorCode result;
 char line[10000];
 file.getline( line, 10000 );
 std::string a = line;
 if( std::string::npos == a.find( "Tally bin boundaries:" ) ) return MB_FAILURE;
 
 // Decide what coordinate system the tally is using
 // first check for Cylindrical coordinates:
 file.getline( line, 10000 );
 a = line;
 std::string::size_type b = a.find( "Cylinder origin at" );
 if( std::string::npos != b )
 {
 coord_sys = CYLINDRICAL;
 if( debug ) std::cout << "origin, axis, direction=| " << a << std::endl;
 std::istringstream ss( a.substr( b + sizeof( "Cylinder origin at" ), 10000 ) );
 // The meshtal file does not contain sufficient information to transform
 // the particle. Although origin, axs, and vec is needed only origin and
 // axs appear in the meshtal file. Why was vec omitted?.
 // get origin (not used)
 // Cylinder origin at 0.00E+00 0.00E+00 0.00E+00,
 // axis in 0.000E+00 0.000E+00 1.000E+00 direction
 double origin[3];
 if( debug ) std::cout << "origin=| ";
 for( int i = 0; i < 3; i++ )
 {
 ss >> origin[i];
 if( debug ) std::cout << origin[i] << " ";
 }
 if( debug ) std::cout << std::endl;
 int length_of_string = 10;
 ss.ignore( length_of_string, ' ' );
 ss.ignore( length_of_string, ' ' );
 ss.ignore( length_of_string, ' ' );
 // Get axis (not used)
 double axis[3];
 if( debug ) std::cout << "axis=| ";
 for( int i = 0; i < 3; i++ )
 {
 ss >> axis[i];
 if( debug ) std::cout << axis[i] << " ";
 }
 if( debug ) std::cout << std::endl;
 file.getline( line, 10000 );
 a = line;
 
 // Get r planes
 if( debug ) std::cout << "R direction:=| ";
 b = a.find( "R direction:" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "R direction" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[0] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Get z planes
 file.getline( line, 10000 );
 a = line;
 if( debug ) std::cout << "Z direction:=| ";
 b = a.find( "Z direction:" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[1] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Get theta planes
 file.getline( line, 10000 );
 a = line;
 if( debug ) std::cout << "Theta direction:=| ";
 b = a.find( "Theta direction (revolutions):" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "Theta direction (revolutions):" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[2] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Cartesian coordinate system:
 }
 else if( std::string::npos != a.find( "X direction:" ) )
 {
 coord_sys = CARTESIAN;
 // Get x planes
 if( debug ) std::cout << "X direction:=| ";
 b = a.find( "X direction:" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "X direction" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[0] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Get y planes
 file.getline( line, 10000 );
 a = line;
 if( debug ) std::cout << "Y direction:=| ";
 b = a.find( "Y direction:" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "Y direction" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[1] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Get z planes
 file.getline( line, 10000 );
 a = line;
 if( debug ) std::cout << "Z direction:=| ";
 b = a.find( "Z direction:" );
 if( std::string::npos != b )
 {
 std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
 result = get_mesh_plane( ss2, debug, planes[2] );
 if( MB_SUCCESS != result ) return result;
 }
 else
 return MB_FAILURE;
 
 // Spherical coordinate system not yet implemented:
 }
 else
 return MB_FAILURE;
 
 return MB_SUCCESS;
}

References CARTESIAN, CYLINDRICAL, moab::debug, ErrorCode, get_mesh_plane(), and MB_SUCCESS.

Referenced by load_one_file().

read_tag_values()

ErrorCode moab::ReadMCNP5::read_tag_values ( const char *  file_name, const char *  tag_name, const FileOptions &  opts, std::vector< int > &  tag_values_out, const SubsetList *  subset_list = 0  )

virtual

Read tag values from a file.

Read the list if all integer tag values from the file for a tag that is a single integer value per entity.

Parameters

file_name	The file to read.
tag_name	The tag for which to read values
tag_values_out	Output: The list of tag values.
subset_list	An array of tag name and value sets specifying the subset of the file to read. If multiple tags are specified, the sets that match all tags (intersection) should be read.
subset_list_length	The length of the 'subset_list' array.

Implements moab::ReaderIface.

Definition at line 63 of file ReadMCNP5.cpp.

{
 return MB_NOT_IMPLEMENTED;
}

References MB_NOT_IMPLEMENTED.

read_tally_header()

ErrorCode moab::ReadMCNP5::read_tally_header ( std::fstream &  file, bool  debug, unsigned int &  tally_number, char  tally_comment[100], particle &  tally_particle  )

private

Definition at line 430 of file ReadMCNP5.cpp.

{
 // Get tally number
 // Mesh Tally Number 104
 ErrorCode result;
 char line[100];
 file.getline( line, 100 );
 std::string a = line;
 std::string::size_type b = a.find( "Mesh Tally Number" );
 if( std::string::npos != b )
 {
 std::istringstream tally_number_ss( a.substr( b + sizeof( "Mesh Tally Number" ), 100 ) );
 tally_number_ss >> tally_number;
 if( debug ) std::cout << "tally_number=| " << tally_number << std::endl;
 }
 else
 {
 std::cout << "tally number not found" << std::endl;
 return MB_FAILURE;
 }
 
 // Next get the tally comment (optional) and particle type
 // 3mm neutron heating in Be (W/cc)
 // This is a neutron mesh tally.
 // std::string tally_comment;
 
 // Get tally particle
 file.getline( line, 100 );
 a = line;
 result = get_tally_particle( a, debug, tally_particle );
 if( MB_FAILURE == result )
 {
 // If this line does not specify the particle type, then it is a tally comment.
 // Get the comment, then get the particle type from the next line.
 tally_comment = line;
 file.getline( line, 100 );
 a = line;
 result = get_tally_particle( a, debug, tally_particle );
 if( MB_SUCCESS != result ) return result;
 }
 if( debug ) std::cout << "tally_comment=| " << tally_comment << std::endl;
 
 return MB_SUCCESS;
}

References moab::debug, ErrorCode, get_tally_particle(), and MB_SUCCESS.

Referenced by load_one_file().

set_header_tags()

ErrorCode moab::ReadMCNP5::set_header_tags ( EntityHandle  output_meshset, char  date_and_time[100], char  title[100], unsigned long int  nps, Tag  data_and_time_tag, Tag  title_tag, Tag  nps_tag  )

private

Definition at line 411 of file ReadMCNP5.cpp.

{
 ErrorCode result;
 result = MBI->tag_set_data( data_and_time_tag, &output_meshset, 1, &date_and_time );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( title_tag, &output_meshset, 1, &title );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( nps_tag, &output_meshset, 1, &nps );
 if( MB_SUCCESS != result ) return result;
 
 return MB_SUCCESS;
}

References ErrorCode, MB_SUCCESS, MBI, and moab::Interface::tag_set_data().

Referenced by load_one_file().

set_tally_tags()

ErrorCode moab::ReadMCNP5::set_tally_tags ( EntityHandle  tally_meshset, unsigned int  tally_number, char  tally_comment[100], particle  tally_particle, coordinate_system  tally_coord_sys, Tag  tally_number_tag, Tag  tally_comment_tag, Tag  tally_particle_tag, Tag  tally_coord_sys_tag  )

private

Definition at line 718 of file ReadMCNP5.cpp.

{
 ErrorCode result;
 result = MBI->tag_set_data( tally_number_tag, &tally_meshset, 1, &tally_number );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( tally_comment_tag, &tally_meshset, 1, &tally_comment );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( tally_particle_tag, &tally_meshset, 1, &tally_particle );
 if( MB_SUCCESS != result ) return result;
 result = MBI->tag_set_data( tally_coord_sys_tag, &tally_meshset, 1, &tally_coord_sys );
 if( MB_SUCCESS != result ) return result;
 
 return MB_SUCCESS;
}

References ErrorCode, MB_SUCCESS, MBI, and moab::Interface::tag_set_data().

Referenced by load_one_file().

transform_point_to_cartesian()

ErrorCode moab::ReadMCNP5::transform_point_to_cartesian ( double *  in, double *  out, coordinate_system  coord_sys  )

private

Definition at line 979 of file ReadMCNP5.cpp.

{
 // Transform coordinate system
 switch( coord_sys )
 {
 case CARTESIAN:
 out[0] = in[0]; // x
 out[1] = in[1]; // y
 out[2] = in[2]; // z
 break;
 // theta is in rotations
 case CYLINDRICAL:
 out[0] = in[0] * cos( 2 * PI * in[2] ); // x
 out[1] = in[0] * sin( 2 * PI * in[2] ); // y
 out[2] = in[1]; // z
 break;
 case SPHERICAL:
 return MB_NOT_IMPLEMENTED;
 default:
 return MB_NOT_IMPLEMENTED;
 }
 
 return MB_SUCCESS;
}

References CARTESIAN, CYLINDRICAL, MB_NOT_IMPLEMENTED, MB_SUCCESS, PI, and SPHERICAL.

Referenced by create_vertices().

Member Data Documentation

C2PI

const double moab::ReadMCNP5::C2PI = 0.1591549430918954

staticprivate

Definition at line 81 of file ReadMCNP5.hpp.

CPI

const double moab::ReadMCNP5::CPI = 0.3183098861837907

staticprivate

Definition at line 82 of file ReadMCNP5.hpp.

elemId

int moab::ReadMCNP5::elemId

private

Definition at line 105 of file ReadMCNP5.hpp.

Referenced by create_elements(), and load_file().

fileIDTag

const Tag* moab::ReadMCNP5::fileIDTag

private

Definition at line 104 of file ReadMCNP5.hpp.

Referenced by create_elements(), create_vertices(), and load_file().

MBI

Interface* moab::ReadMCNP5::MBI

private

Definition at line 102 of file ReadMCNP5.hpp.

Referenced by average_with_existing_tally(), create_elements(), create_tags(), create_vertices(), load_one_file(), ReadMCNP5(), set_header_tags(), set_tally_tags(), and ~ReadMCNP5().

nodeId

int moab::ReadMCNP5::nodeId

private

Definition at line 105 of file ReadMCNP5.hpp.

Referenced by create_vertices(), and load_file().

PI

const double moab::ReadMCNP5::PI = 3.141592653589793

staticprivate

Definition at line 80 of file ReadMCNP5.hpp.

Referenced by transform_point_to_cartesian().

readMeshIface

ReadUtilIface* moab::ReadMCNP5::readMeshIface

private

Definition at line 99 of file ReadMCNP5.hpp.

Referenced by create_elements(), create_vertices(), ReadMCNP5(), and ~ReadMCNP5().

---

The documentation for this class was generated from the following files:

• ReadMCNP5.hpp
• ReadMCNP5.cpp