ReadMCNP5.cpp

Go to the documentation of this file.

#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 <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>
#include "assert.h"
#include "math.h"

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

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

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

// 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( !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 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;

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

// Average a tally that was recently read in with one that alread 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) 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 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(!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