moab
|
#include <ReadMCNP5.hpp>
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 51 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 88 of file ReadMCNP5.hpp.
moab::ReadMCNP5::ReadMCNP5 | ( | Interface * | impl = NULL | ) |
Definition at line 43 of file ReadMCNP5.cpp.
: MBI(impl), fileIDTag(0) { assert( NULL!=impl); MBI->query_interface(readMeshIface); assert( NULL!=readMeshIface ); }
moab::ReadMCNP5::~ReadMCNP5 | ( | ) | [virtual] |
Definition at line 51 of file ReadMCNP5.cpp.
{ 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 983 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(!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 882 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) return result; result = MBI->tag_get_data( error_tag, existing_elements, errors0 ); if(MB_SUCCESS != result) return result; // average the values and errors result = average_tally_values( nps0, nps1, values0, values1, errors0, errors1, n_elements ); if(MB_SUCCESS != result) 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) return result; result = MBI->tag_set_data( error_tag, existing_elements, errors0 ); if(MB_SUCCESS != result) 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 797 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 connectivitiy 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 389 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; }
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 748 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; 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 38 of file ReadMCNP5.cpp.
ErrorCode moab::ReadMCNP5::get_mesh_plane | ( | std::istringstream & | ss, |
bool | debug, | ||
std::vector< double > & | plane | ||
) | [private] |
Definition at line 662 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; }
ErrorCode moab::ReadMCNP5::get_tally_particle | ( | std::string | a, |
bool | debug, | ||
particle & | tally_particle | ||
) | [private] |
Definition at line 524 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; }
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 70 of file ReadMCNP5.cpp.
{ // at this time there is no support for reading a subset of the file if (subset_list) { readMeshIface->report_error( "Reading subset of files not supported for meshtal." ); return MB_UNSUPPORTED_OPERATION; } nodeId = elemId = 0; fileIDTag = file_id_tag; // Average several meshtal files if the AVERAGE_TALLY option is givin. // 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 133 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); // 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 cylidrical 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 cylidrical 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 676 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( !finite(errors[index]) ) { std::cerr << "found nan error while reading file" << std::endl; errors[index] = 1.0; } if( !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 429 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; }
ErrorCode moab::ReadMCNP5::read_mesh_planes | ( | std::fstream & | file, |
bool | debug, | ||
std::vector< double > | planes[3], | ||
coordinate_system & | coord_sys | ||
) | [private] |
Definition at line 540 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; }
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 59 of file ReadMCNP5.cpp.
{ 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 480 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; }
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 463 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; }
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 727 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; }
ErrorCode moab::ReadMCNP5::transform_point_to_cartesian | ( | double * | in, |
double * | out, | ||
coordinate_system | coord_sys | ||
) | [private] |
Definition at line 957 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; }
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 99 of file ReadMCNP5.hpp.
const Tag* moab::ReadMCNP5::fileIDTag [private] |
Definition at line 98 of file ReadMCNP5.hpp.
Interface* moab::ReadMCNP5::MBI [private] |
Definition at line 96 of file ReadMCNP5.hpp.
int moab::ReadMCNP5::nodeId [private] |
Definition at line 99 of file ReadMCNP5.hpp.
const double moab::ReadMCNP5::PI = 3.141592653589793 [static, private] |
Definition at line 80 of file ReadMCNP5.hpp.
ReadUtilIface* moab::ReadMCNP5::readMeshIface [private] |
Definition at line 93 of file ReadMCNP5.hpp.