MOAB: Mesh Oriented datABase  (version 5.4.1)
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.
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 ReaderIfacefactory (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

ReadUtilIfacereadMeshIface
InterfaceMBI
const TagfileIDTag
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

Enumerator:
NO_SYSTEM 
CARTESIAN 
CYLINDRICAL 
SPHERICAL 

Definition at line 84 of file ReadMCNP5.hpp.

enum moab::ReadMCNP5::particle [private]
Enumerator:
NEUTRON 
PHOTON 
ELECTRON 

Definition at line 91 of file ReadMCNP5.hpp.


Constructor & Destructor Documentation

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 );
}

Member Function Documentation

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;
}

Definition at line 40 of file ReadMCNP5.cpp.

References ReadMCNP5().

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

{
    return new ReadMCNP5( iface );
}
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.

Parameters:
file_nameThe file to read.
file_setOptional 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_listAn optional struct pointer specifying the tags identifying entity sets to be read.
file_id_tagIf 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.

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(), MB_SUCCESS, and origin.

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.

Parameters:
file_nameThe file to read.
tag_nameThe tag for which to read values
tag_values_outOutput: The list of tag values.
subset_listAn 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_lengthThe 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;
}

Member Data Documentation

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().

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().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines