 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#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. | |
| 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. | |
| ReadMCNP5 (Interface *impl=NULL) | |
| virtual | ~ReadMCNP5 () |
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 |
Definition at line 52 of file ReadMCNP5.hpp.
enum moab::ReadMCNP5::coordinate_system [private] |
Definition at line 84 of file ReadMCNP5.hpp.
{
NO_SYSTEM,
CARTESIAN,
CYLINDRICAL,
SPHERICAL
};
enum moab::ReadMCNP5::particle [private] |
Definition at line 91 of file ReadMCNP5.hpp.
| moab::ReadMCNP5::ReadMCNP5 | ( | Interface * | impl = NULL | ) |
Definition at line 46 of file ReadMCNP5.cpp.
References MBI, moab::Interface::query_interface(), and readMeshIface.
Referenced by factory().
: MBI( impl ), fileIDTag( NULL ), nodeId( 0 ), elemId( 0 ) { assert( NULL != impl ); MBI->query_interface( readMeshIface ); assert( NULL != readMeshIface ); }
| moab::ReadMCNP5::~ReadMCNP5 | ( | ) | [virtual] |
Definition at line 54 of file ReadMCNP5.cpp.
References MBI, readMeshIface, and moab::Interface::release_interface().
{
if( readMeshIface )
{
MBI->release_interface( readMeshIface );
readMeshIface = 0;
}
}
| 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.
References moab::Util::is_finite(), and MB_SUCCESS.
Referenced by average_with_existing_tally().
{
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;
}
| 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.
References average_tally_values(), 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().
{
// 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;
}
| 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.
References moab::Interface::add_entities(), moab::ReadUtilIface::assign_ids(), CARTESIAN, CYLINDRICAL, 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().
{
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;
}
| 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.
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().
{
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;
}
| 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.
References moab::Interface::add_entities(), moab::ReadUtilIface::assign_ids(), 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().
{
// 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;
}
| ReaderIface * moab::ReadMCNP5::factory | ( | Interface * | iface | ) | [static] |
Definition at line 40 of file ReadMCNP5.cpp.
References ReadMCNP5().
Referenced by moab::ReaderWriterSet::ReaderWriterSet().
| ErrorCode moab::ReadMCNP5::get_mesh_plane | ( | std::istringstream & | ss, |
| bool | debug, | ||
| std::vector< double > & | plane | ||
| ) | [private] |
Definition at line 647 of file ReadMCNP5.cpp.
References MB_SUCCESS.
Referenced by read_mesh_planes().
{
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;
}
| ErrorCode moab::ReadMCNP5::get_tally_particle | ( | std::string | a, |
| bool | debug, | ||
| particle & | tally_particle | ||
| ) | [private] |
Definition at line 479 of file ReadMCNP5.cpp.
References ELECTRON, MB_SUCCESS, NEUTRON, and PHOTON.
Referenced by read_tally_header().
{
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;
}
| 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.
| 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. |
Implements moab::ReaderIface.
Definition at line 73 of file ReadMCNP5.cpp.
References elemId, ErrorCode, fileIDTag, moab::FileOptions::get_int_option(), length(), load_one_file(), MB_SET_ERR, MB_SUCCESS, MB_UNSUPPORTED_OPERATION, and nodeId.
{
// 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;
}
| 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.
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().
{
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;
}
| 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.
References moab::Util::is_finite(), MB_SUCCESS, PHOTON, and size.
Referenced by load_one_file().
{
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;
}
| 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.
References MB_SUCCESS.
Referenced by load_one_file().
{
// 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;
}
| 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.
References CARTESIAN, CYLINDRICAL, ErrorCode, get_mesh_plane(), and MB_SUCCESS.
Referenced by load_one_file().
{
// 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;
}
| 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.
| 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.
References MB_NOT_IMPLEMENTED.
{
return MB_NOT_IMPLEMENTED;
}
| 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.
References ErrorCode, get_tally_particle(), and MB_SUCCESS.
Referenced by load_one_file().
{
// 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;
}
| 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.
References ErrorCode, MB_SUCCESS, MBI, and moab::Interface::tag_set_data().
Referenced by load_one_file().
{
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;
}
| 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.
References ErrorCode, MB_SUCCESS, MBI, and moab::Interface::tag_set_data().
Referenced by load_one_file().
{
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;
}
| ErrorCode moab::ReadMCNP5::transform_point_to_cartesian | ( | double * | in, |
| double * | out, | ||
| coordinate_system | coord_sys | ||
| ) | [private] |
Definition at line 979 of file ReadMCNP5.cpp.
References CARTESIAN, CYLINDRICAL, MB_NOT_IMPLEMENTED, MB_SUCCESS, PI, and SPHERICAL.
Referenced by create_vertices().
{
// 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;
}
const double moab::ReadMCNP5::C2PI = 0.1591549430918954 [static, private] |
Definition at line 81 of file ReadMCNP5.hpp.
const double moab::ReadMCNP5::CPI = 0.3183098861837907 [static, private] |
Definition at line 82 of file ReadMCNP5.hpp.
int moab::ReadMCNP5::elemId [private] |
Definition at line 105 of file ReadMCNP5.hpp.
Referenced by create_elements(), and load_file().
const Tag* moab::ReadMCNP5::fileIDTag [private] |
Definition at line 104 of file ReadMCNP5.hpp.
Referenced by create_elements(), create_vertices(), and load_file().
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().
int moab::ReadMCNP5::nodeId [private] |
Definition at line 105 of file ReadMCNP5.hpp.
Referenced by create_vertices(), and load_file().
const double moab::ReadMCNP5::PI = 3.141592653589793 [static, private] |
Definition at line 80 of file ReadMCNP5.hpp.
Referenced by transform_point_to_cartesian().
ReadUtilIface* moab::ReadMCNP5::readMeshIface [private] |
Definition at line 99 of file ReadMCNP5.hpp.
Referenced by create_elements(), create_vertices(), ReadMCNP5(), and ~ReadMCNP5().
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.