ReadMCNP5.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
#include "ReadMCNP5.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "Internals.hpp"  // For MB_START_ID
#include "moab/Range.hpp"
#include "moab/FileOptions.hpp"
#include "moab/Util.hpp"
 
#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <cassert>
 
namespace moab
{
 
// Parameters
const double ReadMCNP5::PI   = 3.141592653589793;
const double ReadMCNP5::C2PI = 0.1591549430918954;
const double ReadMCNP5::CPI  = 0.3183098861837907;
 
ReaderIface* ReadMCNP5::factory( Interface* iface )
{
    return new ReadMCNP5( iface );
}
 
// Constructor
ReadMCNP5::ReadMCNP5( Interface* impl ) : MBI( impl ), fileIDTag( NULL ), nodeId( 0 ), elemId( 0 )
{
    assert( NULL != impl );
    MBI->query_interface( readMeshIface );
    assert( NULL != readMeshIface );
}
 
// Destructor
ReadMCNP5::~ReadMCNP5()
{
    if( readMeshIface )
    {
        MBI->release_interface( readMeshIface );
        readMeshIface = 0;
    }
}
 
ErrorCode 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 */ )
{
    return MB_NOT_IMPLEMENTED;
}
 
// Load the file as called by the Interface function
ErrorCode ReadMCNP5::load_file( const char* filename,
                                const EntityHandle* input_meshset,
                                const FileOptions& options,
                                const ReaderIface::SubsetList* subset_list,
                                const Tag* file_id_tag )
{
    // 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;
}
 
// This actually reads the file. It creates the mesh elements unless
// the file is being averaged with a pre-existing mesh.
ErrorCode ReadMCNP5::load_one_file( const char* fname,
                                    const EntityHandle* input_meshset,
                                    const FileOptions& options,
                                    const bool average )
{
    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;
}
 
// create tags needed for this reader
ErrorCode 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 )
{
    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 ReadMCNP5::read_file_header( std::fstream& file,
                                       bool debug,
                                       char date_and_time[100],
                                       char title[100],
                                       unsigned long int& nps )
{
    // 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 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 )
{
    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 ReadMCNP5::read_tally_header( std::fstream& file,
                                        bool debug,
                                        unsigned int& tally_number,
                                        char tally_comment[100],
                                        particle& tally_particle )
{
    // 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 ReadMCNP5::get_tally_particle( std::string a, bool debug, particle& tally_particle )
{
    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 ReadMCNP5::read_mesh_planes( std::fstream& file,
                                       bool debug,
                                       std::vector< double > planes[3],
                                       coordinate_system& coord_sys )
{
    // 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;
}
 
// Given a stringstream, return a vector of values in the string.
ErrorCode ReadMCNP5::get_mesh_plane( std::istringstream& ss, bool debug, std::vector< double >& plane )
{
    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 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[] )
{
    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 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 )
{
    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 ReadMCNP5::create_vertices( std::vector< double > planes[3],
                                      bool debug,
                                      EntityHandle& start_vert,
                                      coordinate_system coord_sys,
                                      EntityHandle tally_meshset )
{
    // 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;
}
 
ErrorCode 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 )
{
    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;
}
 
// Average a tally that was recently read in with one that already exists in
// the interface. Only the existing values will be updated.
ErrorCode ReadMCNP5::average_with_existing_tally( bool debug,
                                                  unsigned long int& new_nps,
                                                  unsigned long int nps1,
                                                  unsigned int tally_number,
                                                  Tag tally_number_tag,
                                                  Tag nps_tag,
                                                  Tag tally_tag,
                                                  Tag error_tag,
                                                  double* values1,
                                                  double* errors1,
                                                  unsigned int n_elements )
{
    // 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 ReadMCNP5::transform_point_to_cartesian( double* in, double* out, coordinate_system coord_sys )
{
    // 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;
}
 
// Average two tally values and their error. Return average values in the
// place of first tally values.
ErrorCode 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 )
{
    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;
}
 
}  // namespace moab