moab
moab::ReadMCNP5 Class Reference

#include <ReadMCNP5.hpp>

Inheritance diagram for moab::ReadMCNP5:
moab::ReaderIface

List of all members.

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 51 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 88 of file ReadMCNP5.hpp.


Constructor & Destructor Documentation

Definition at line 43 of file ReadMCNP5.cpp.

  : MBI(impl), fileIDTag(0) {
    assert( NULL!=impl);
    MBI->query_interface(readMeshIface);
    assert( NULL!=readMeshIface );
}

Definition at line 51 of file ReadMCNP5.cpp.


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 983 of file ReadMCNP5.cpp.

                                                                              {
  
  for(unsigned long int i=0; i<n_values; i++) {
    //std::cout << " values0=" << values0[i] << " values1=" << values1[i] 
    //          << " errors0=" << errors0[i] << " errors1=" << errors1[i] 
    //          << " nps0=" << nps0 << " nps1=" << nps1 << std::endl;
    errors0[i] = sqrt( pow(values0[i]*errors0[i]*nps0,2) + 
                       pow(values1[i]*errors1[i]*nps1,2) ) / 
                 (values0[i]*nps0 + values1[i]*nps1);

    // It is possible to introduce nans if the values are zero.
    if(!finite(errors0[i])) errors0[i] = 1.0;

    values0[i] = ( values0[i]*nps0 + values1[i]*nps1 ) / (nps0 + nps1);

    //std::cout << " values0=" << values0[i] << " errors0=" << errors0[i] << std::endl;
  }
  // REMEMBER TO UPDATE NPS0 = NPS0 + NPS1 after this
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::average_with_existing_tally ( bool  debug,
unsigned long int &  new_nps,
unsigned long int  nps,
unsigned int  tally_number,
Tag  tally_number_tag,
Tag  nps_tag,
Tag  tally_tag,
Tag  error_tag,
double  values[],
double  errors[],
unsigned int  n_elements 
) [private]

Definition at line 882 of file ReadMCNP5.cpp.

                                                                              {
    
  // get the tally number
  ErrorCode result;

  // match the tally number with one from the existing meshtal file
  Range matching_tally_number_sets;
  const void* const tally_number_val[] = {&tally_number};
  result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &tally_number_tag, 
                                              tally_number_val, 1, 
                                              matching_tally_number_sets );
  if(MB_SUCCESS != result) return result;
  if (debug) std::cout << "number of matching meshsets=" 
                       << matching_tally_number_sets.size() << std::endl;
  assert(1 == matching_tally_number_sets.size());

  // identify which of the meshsets is existing
  EntityHandle existing_meshset;
  existing_meshset = matching_tally_number_sets.front();

  // get the existing elements from the set
  Range existing_elements;
  result = MBI->get_entities_by_type( existing_meshset, MBHEX, existing_elements );
  if(MB_SUCCESS != result) return result;
  assert(existing_elements.size() == n_elements);

  // get the nps of the existing and new tally
  unsigned long int nps0;
  Range sets_with_this_tag;
  result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, 
                                              sets_with_this_tag);
  if(MB_SUCCESS != result) return result;
  if (debug) std::cout << "number of nps sets=" << sets_with_this_tag.size() << std::endl;
  assert(1 == sets_with_this_tag.size());
  result = MBI->tag_get_data( nps_tag, &sets_with_this_tag.front(), 1, &nps0);
  if(MB_SUCCESS != result) return result;
  if (debug) std::cout << "nps0=" << nps0 << " nps1=" << nps1 << std::endl;
  new_nps = nps0 + nps1;

  // get tally values from the existing elements
  double *values0 = new double [existing_elements.size()];
  double *errors0 = new double [existing_elements.size()];
  result = MBI->tag_get_data( tally_tag, existing_elements, values0 );
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_get_data( error_tag, existing_elements, errors0 );
  if(MB_SUCCESS != result) return result;

  // average the values and errors
  result = average_tally_values( nps0, nps1, values0, values1, 
                                 errors0, errors1, n_elements );
  if(MB_SUCCESS != result) return result;
  
  // set the averaged information back onto the existing elements
  result = MBI->tag_set_data( tally_tag, existing_elements, values0 );
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( error_tag, existing_elements, errors0 );
  if(MB_SUCCESS != result) return result;
  
  // cleanup
  delete[] values0;
  delete[] errors0;

  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::create_elements ( bool  debug,
std::vector< double >  planes[3],
unsigned int  n_chopped_x0_planes,
unsigned int  n_chopped_x2_planes,
EntityHandle  start_vert,
double  values[],
double  errors[],
Tag  tally_tag,
Tag  error_tag,
EntityHandle  tally_meshset,
coordinate_system  tally_coord_sys 
) [private]

Definition at line 797 of file ReadMCNP5.cpp.

                                                                              {
  ErrorCode result;
  unsigned int index;
  EntityHandle start_element = 0;
  unsigned int n_elements = (planes[0].size()-1) * (planes[1].size()-1) 
                                                 * (planes[2].size()-1);
  EntityHandle *connect;
  result = readMeshIface->get_element_connect( n_elements, 8, MBHEX, MB_START_ID, 
                                             start_element, connect );
  if(MB_SUCCESS != result) return result;
  assert(0 != start_element); // check for NULL

  unsigned int counter = 0;
  for (unsigned int i=0; i<planes[0].size()-1; i++) {
    for (unsigned int j=0; j<planes[1].size()-1; j++) {
      for (unsigned int k=0; k<planes[2].size()-1; k++) {
          
        index = start_vert + i + j*planes[0].size() + k*planes[0].size()*planes[1].size();
        // for rectangular mesh, the file prints: x y z
        // z changes the fastest and x changes the slowest.
        // This means that the connectivitiy ordering is not consistent between
        // rectangular and cylindrical mesh.
        if(CARTESIAN == tally_coord_sys) {
          connect[0] = index;
          connect[1] = index + 1;  
          connect[2] = index + 1 + planes[0].size();     
          connect[3] = index +     planes[0].size();
          connect[4] = index +                        planes[0].size()*planes[1].size();
          connect[5] = index + 1 +                    planes[0].size()*planes[1].size();    
          connect[6] = index + 1 + planes[0].size() + planes[0].size()*planes[1].size();
          connect[7] = index +     planes[0].size() + planes[0].size()*planes[1].size();        
        // for cylindrical mesh, the file prints: r z theta
        // Theta changes the fastest and r changes the slowest.
        } else if(CYLINDRICAL == tally_coord_sys) {
          connect[0] = index;
          connect[1] = index + 1;  
          connect[2] = index + 1 +                    planes[0].size()*planes[1].size();     
          connect[3] = index +                        planes[0].size()*planes[1].size();
          connect[4] = index +     planes[0].size();
          connect[5] = index + 1 + planes[0].size();    
          connect[6] = index + 1 + planes[0].size() + planes[0].size()*planes[1].size();
          connect[7] = index +     planes[0].size() + planes[0].size()*planes[1].size();
    } else return MB_NOT_IMPLEMENTED;

    connect += 8;
    counter++;
      }
    }
  }
  if (counter != n_elements) std::cout << "counter=" << counter << " n_elements=" 
                                       << n_elements << std::endl;

  Range element_range(start_element, start_element+n_elements-1);
  result = MBI->tag_set_data(tally_tag, element_range, values);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data(error_tag, element_range, errors);
  if(MB_SUCCESS != result) return result;
  
  // add the elements to the tally set
  result = MBI->add_entities( tally_meshset, element_range );
  if(MB_SUCCESS != result) return result;
  if (debug) std::cout << "Read " << n_elements << " elements from tally." << std::endl;
  
  if (fileIDTag) {
    result = readMeshIface->assign_ids( *fileIDTag, element_range, elemId );
    if (MB_SUCCESS != result)
      return result;
    elemId += element_range.size();
  }

  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::create_tags ( Tag date_and_time_tag,
Tag title_tag,
Tag nps_tag,
Tag tally_number_tag,
Tag tally_comment_tag,
Tag tally_particle_tag,
Tag tally_coord_sys_tag,
Tag tally_tag,
Tag error_tag 
) [private]

Definition at line 389 of file ReadMCNP5.cpp.

                                                     {
  ErrorCode result;
  result = MBI->tag_get_handle("DATE_AND_TIME_TAG", 100, MB_TYPE_OPAQUE, 
                           date_and_time_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TITLE_TAG", 100, MB_TYPE_OPAQUE, title_tag,
                           MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("NPS_TAG", sizeof(unsigned long int), MB_TYPE_OPAQUE,
                           nps_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TALLY_NUMBER_TAG", 1, MB_TYPE_INTEGER, 
                           tally_number_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TALLY_COMMENT_TAG", 100, MB_TYPE_OPAQUE, 
                           tally_comment_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TALLY_PARTICLE_TAG", sizeof(particle), MB_TYPE_OPAQUE,
                           tally_particle_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TALLY_COORD_SYS_TAG", sizeof(coordinate_system), MB_TYPE_OPAQUE,
                           tally_coord_sys_tag, MB_TAG_SPARSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("TALLY_TAG", 1, MB_TYPE_DOUBLE, tally_tag,
                           MB_TAG_DENSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  result = MBI->tag_get_handle("ERROR_TAG", 1, MB_TYPE_DOUBLE, error_tag,
                           MB_TAG_DENSE|MB_TAG_CREAT);
  if(MB_SUCCESS!=result) return result;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::create_vertices ( std::vector< double >  planes[3],
bool  debug,
EntityHandle start_vert,
coordinate_system  coord_sys,
EntityHandle  tally_meshset 
) [private]

Definition at line 748 of file ReadMCNP5.cpp.

                                                                          {
                                         
  // The only info needed to build elements is the mesh plane boundaries.
  ErrorCode result;
  int n_verts = planes[0].size() * planes[1].size() * planes[2].size();
  if (debug) std::cout << "n_verts=" << n_verts << std::endl;
  std::vector<double*> coord_arrays(3);
  result = readMeshIface->get_node_coords( 3, n_verts, MB_START_ID, 
                                           start_vert, coord_arrays );
  if(MB_SUCCESS != result) return result;
  assert(0 != start_vert); // check for NULL

  for (unsigned int k=0; k < planes[2].size(); k++) {
    for (unsigned int j=0; j < planes[1].size(); j++) {
      for (unsigned int i=0; i < planes[0].size(); i++) {

        unsigned int idx = (k*planes[0].size()*planes[1].size() + j*planes[0].size() + i);
        double in[3], out[3];

        in[0] = planes[0][i];
        in[1] = planes[1][j];
        in[2] = planes[2][k];
        result = transform_point_to_cartesian( in, out, coord_sys );
        if(MB_SUCCESS != result) return result;

        coord_arrays[0][idx] = out[0];
        coord_arrays[1][idx] = out[1];
        coord_arrays[2][idx] = out[2];
      }
    }
  }
  Range vert_range(start_vert, start_vert+n_verts-1);
  result = MBI->add_entities( tally_meshset, vert_range );
  if(MB_SUCCESS != result) return result;
  
  if (fileIDTag) {
    result = readMeshIface->assign_ids( *fileIDTag, vert_range, nodeId );
    if (MB_SUCCESS != result)
      return result;
    nodeId += vert_range.size();
  }
  
  return MB_SUCCESS;
}

Definition at line 38 of file ReadMCNP5.cpp.

                                                  { 
  return new ReadMCNP5( iface );
}
ErrorCode moab::ReadMCNP5::get_mesh_plane ( std::istringstream &  ss,
bool  debug,
std::vector< double > &  plane 
) [private]

Definition at line 662 of file ReadMCNP5.cpp.

                                                                 {
  double value;
  plane.clear();
  while (!ss.eof()) {
    ss >> value;
    plane.push_back( value );
    if (debug) std::cout << value << " ";
  }
  if (debug) std::cout << std::endl;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::get_tally_particle ( std::string  a,
bool  debug,
particle tally_particle 
) [private]

Definition at line 524 of file ReadMCNP5.cpp.

                                                                            {
  
  if        (std::string::npos != a.find("This is a neutron mesh tally.")) {
    tally_particle = NEUTRON;
  } else if (std::string::npos != a.find("This is a photon mesh tally.")) {
    tally_particle = PHOTON;
  } else if (std::string::npos != a.find("This is an electron mesh tally.")) {
    tally_particle = ELECTRON;
  } else return MB_FAILURE;

  if (debug) std::cout << "tally_particle=| " << tally_particle << std::endl;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::load_file ( const char *  file_name,
const EntityHandle file_set,
const FileOptions opts,
const SubsetList subset_list = 0,
const Tag file_id_tag = 0 
) [virtual]

Load mesh from a file.

Method all readers must provide to import a mesh.

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 70 of file ReadMCNP5.cpp.

                                                                           {
  // at this time there is no support for reading a subset of the file
  if (subset_list) {
    readMeshIface->report_error( "Reading subset of files not supported for meshtal." );
    return MB_UNSUPPORTED_OPERATION;
  }
  
  nodeId = elemId = 0;
  fileIDTag = file_id_tag;

  // Average several meshtal files if the AVERAGE_TALLY option is givin.
  // In this case, the integer value is the number of files to average.
  // If averaging, the filename passed to load_file is assumed to be the
  // root filename. The files are indexed as "root_filename""index".meshtal.
  // Indices start with 1.
  int n_files;
  bool average = false;
  ErrorCode result;
  if(MB_SUCCESS == options.get_int_option("AVERAGE_TALLY", n_files)) {
   
    // read the first file (but do not average -> can't average a single file)
    result = load_one_file( filename, 
                            input_meshset, 
                            options, 
                            average );
    if(MB_SUCCESS != result) return result;

    // get the root filename
    std::string root_filename(filename);
    int length = root_filename.length();
    root_filename.erase(length - sizeof(".meshtal"));

    // average the first file with the rest of the files
    average = true;
    for(int i=2; i<=n_files; i++) {
      std::stringstream index;
      index << i;
      std::string subsequent_filename = root_filename + index.str() + ".meshtal";
      result = load_one_file( subsequent_filename.c_str(), 
                              input_meshset, 
                              options, 
                              average );
      if(MB_SUCCESS != result) return result;
    }
 
  // if not averaging, read a single file
  } else {
    result = load_one_file( filename, 
                            input_meshset, 
                            options, 
                            average );
    if(MB_SUCCESS != result) return result;
  }
  
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::load_one_file ( const char *  fname,
const EntityHandle input_meshset,
const FileOptions options,
const bool  average 
) [private]

Definition at line 133 of file ReadMCNP5.cpp.

                                                                    {

  bool debug = false;
  if (debug) std::cout << "begin ReadMCNP5::load_one_file" << std::endl;

  ErrorCode result;
  std::fstream file;
  file.open( fname, std::fstream::in );
  char line[10000];

  // create tags
  Tag date_and_time_tag,  
        title_tag,           
        nps_tag,  
        tally_number_tag,    
        tally_comment_tag, 
        tally_particle_tag, 
        tally_coord_sys_tag, 
        tally_tag, 
        error_tag;

  result = create_tags( date_and_time_tag,   
                        title_tag,         
                        nps_tag,            
                        tally_number_tag,    
                        tally_comment_tag, 
                        tally_particle_tag, 
                        tally_coord_sys_tag, 
                        tally_tag,         
                        error_tag );
  if(MB_SUCCESS != result) return result;

  // ******************************************************************
  // This info exists only at the top of each meshtal file
  // ******************************************************************

  // define characteristics of the entire file
  char date_and_time[100] = "";
  char title[100] = "";
  // this file's number of particles
  unsigned long int nps;
  // sum of this file's and existing file's nps for averaging
  unsigned long int new_nps;

  // read the file header
  result = read_file_header( file, 
                             debug, 
                             date_and_time, 
                             title,
                             nps );
  if(MB_SUCCESS != result) return result;

  // blank line
  file.getline(line, 10000);

  // Everything stored in the file being read will be in the input_meshset.
  // if this is being saved in MOAB, set header tags
  if (!average && 0 != input_meshset) {
    result = set_header_tags( *input_meshset, 
                              date_and_time,
                              title,
                              nps,
                              date_and_time_tag,
                              title_tag,
                              nps_tag );
    if(MB_SUCCESS != result) return result;
  }

  // ******************************************************************
  // This info is repeated for each tally in the meshtal file.
  // ******************************************************************

  // If averaging, nps will hold the sum of particles simulated in both tallies.
  while( !file.eof() ) {

    // define characteristics of this tally
    unsigned int        tally_number;
    char                tally_comment[100] = "";
    particle            tally_particle;
    coordinate_system   tally_coord_sys;
    std::vector<double> planes[3];
    unsigned int        n_chopped_x0_planes; 
    unsigned int        n_chopped_x2_planes;
    
    // read tally header
    result = read_tally_header( file, 
                                debug, 
                                tally_number, 
                                tally_comment,
                                tally_particle );
    if(MB_SUCCESS != result) return result;
    
    // blank line
    file.getline(line, 10000);

    // read mesh planes
    result = read_mesh_planes( file, 
                               debug, 
                               planes, 
                               tally_coord_sys );
    if(MB_SUCCESS != result) return result;

    // get energy boundaries
    file.getline(line, 10000);
    std::string a = line;
    if (debug) std::cout << "Energy bin boundaries:=| " << a << std::endl;

    // blank
    file.getline(line, 10000);

    // column headers
    file.getline(line, 10000);

    // If using cylidrical mesh, it may be necessary to chop off the last theta element.
    // We chop off the last theta plane because the elements will be wrong and skew up
    // the tree building code. This is
    // because the hex elements are a linear approximation to the cylindrical elements.
    // Chopping off the last plane is problem-dependent, and due to MCNP5's mandate 
    // that the cylidrical mesh must span 360 degrees.
    if (CYLINDRICAL==tally_coord_sys &&
        MB_SUCCESS ==options.get_null_option("REMOVE_LAST_AZIMUTHAL_PLANE")) {
      planes[2].pop_back();
      n_chopped_x2_planes = 1;
      if (debug) std::cout << "remove last cylindrical plane option found" << std::endl;
    } else {
      n_chopped_x2_planes = 0;
    }
    
    // If using cylindrical mesh, it may be necessary to chop off the first radial element.
    // These elements extend from the axis and often do not cover areas of interest. For 
    // example, if the mesh was meant to cover r=390-400, the first layer will go from
    // 0-390 and serve as incorrect source elements for interpolation.  
    if (CYLINDRICAL==tally_coord_sys &&   
        MB_SUCCESS ==options.get_null_option("REMOVE_FIRST_RADIAL_PLANE")) {   
      std::vector<double>::iterator front=planes[0].begin();
      planes[0].erase( front );
      n_chopped_x0_planes = 1;   
      if (debug) std::cout << "remove first radial plane option found" << std::endl;  
    } else {          
      n_chopped_x0_planes = 0;  
    }  

    // Read the values and errors of each element from the file.
    // Do not read values that are chopped off.
    unsigned int n_elements = (planes[0].size()-1)*(planes[1].size()-1)*(planes[2].size()-1);
    double *values = new double [n_elements];
    double *errors = new double [n_elements];
    result = read_element_values_and_errors( file, 
                                             debug, 
                                             planes,
                                             n_chopped_x0_planes,
                                             n_chopped_x2_planes,
                                             tally_particle,
                                             values, 
                                             errors );
    if(MB_SUCCESS != result) return result;
    
    // blank line
    file.getline(line, 10000);
   
    // ****************************************************************
    // This tally has been read. If it is not being averaged, build tags,
    // vertices and elements. If it is being averaged, average the data 
    // with a tally already existing in the MOAB instance.
    // ****************************************************************
    if (!average) {
      EntityHandle tally_meshset;
      result = MBI->create_meshset(MESHSET_SET, tally_meshset);
      if(MB_SUCCESS != result) return result;
      
      // set tags on the tally
      result = set_tally_tags( tally_meshset,
                               tally_number,
                               tally_comment,
                               tally_particle,
                               tally_coord_sys,
                               tally_number_tag, 
                               tally_comment_tag,
                               tally_particle_tag,
                               tally_coord_sys_tag );
      if(MB_SUCCESS != result) return result;

      // The only info needed to build elements is the mesh plane boundaries.
      // Build vertices...
      EntityHandle start_vert = 0;
      result = create_vertices( planes, 
                                debug,
                                start_vert, 
                                tally_coord_sys,
                                tally_meshset );
      if(MB_SUCCESS != result) return result; 
      
      // Build elements and tag them with tally values and errors, then add
      // them to the tally_meshset.
      result = create_elements( debug, 
                                planes,
                                n_chopped_x0_planes, 
                                n_chopped_x2_planes,
                                start_vert, 
                                values, 
                                errors, 
                                tally_tag, 
                                error_tag, 
                                tally_meshset,
                                tally_coord_sys );
      if(MB_SUCCESS != result) return result; 
      
      // add this tally's meshset to the output meshset
      if (debug) std::cout << "not averaging tally" << std::endl;

    // average the tally values, then delete stuff that was created
    } else {
      if (debug) std::cout << "averaging tally" << std::endl;
      result = average_with_existing_tally( debug, 
                                            new_nps,
                                            nps,
                                            tally_number,
                                            tally_number_tag,          
                                            nps_tag,
                                            tally_tag,
                                            error_tag,
                                            values,
                                            errors,
                                            n_elements );
      if(MB_SUCCESS != result) return result;
    }

    // clean up
    delete[] values;
    delete[] errors;
  }

  // If we are averaging, delete the remainder of this file's information.
  // Add the new nps to the existing file's nps if we are averaging.
  // This is calculated during every tally averaging but only used after the last one.
  if (average) {
    Range matching_nps_sets;
    result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 
                                                0, 1, matching_nps_sets );
    if(MB_SUCCESS != result) return result;
    if (debug) std::cout << "number of matching nps meshsets=" 
                         << matching_nps_sets.size() << std::endl;
    assert(1 == matching_nps_sets.size());
    result = MBI->tag_set_data( nps_tag, matching_nps_sets, &new_nps );
    if(MB_SUCCESS != result) return result;

  // If this file is not being averaged, return the output meshset.
  } 
  file.close();
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::read_element_values_and_errors ( std::fstream &  file,
bool  debug,
std::vector< double >  planes[3],
unsigned int  n_chopped_x0_planes,
unsigned int  n_chopped_x2_planes,
particle  tally_particle,
double  values[],
double  errors[] 
) [private]

Definition at line 676 of file ReadMCNP5.cpp.

                                                                           {

  unsigned int index = 0;
  // need to read every line in the file, even if we chop off some elements 
  for (unsigned int i=0; i<planes[0].size()-1+n_chopped_x0_planes; i++) {
    for (unsigned int j=0; j<planes[1].size()-1; j++) {
      for (unsigned int k=0; k<planes[2].size()-1+n_chopped_x2_planes; k++) {
        char line[100];
        file.getline(line, 100);
        // if this element has been chopped off, skip it
        if (i<n_chopped_x0_planes) continue;
        if (k>=planes[2].size()-1 && k<planes[2].size()-1+n_chopped_x2_planes) continue;
        std::string a=line;
        std::stringstream ss(a);
        double centroid[3];
        double energy;
        // for some reason, photon tallies print the energy in the first column
        if (PHOTON==tally_particle) ss >> energy;
        // the centroid is not used in this reader
        ss >> centroid[0];
        ss >> centroid[1];
        ss >> centroid[2];
        // only the tally values and errors are used
        ss >> values[index];
        ss >> errors[index];

        // make sure that input data is good
    if( !finite(errors[index]) ) {
      std::cerr << "found nan error while reading file" << std::endl;
          errors[index] = 1.0;
        }
    if( !finite(values[index]) ) {
      std::cerr << "found nan value while reading file" << std::endl;
          values[index] = 0.0;
        }

        index++;
      }
    }
  }
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::read_file_header ( std::fstream &  file,
bool  debug,
char  date_and_time[100],
char  title[100],
unsigned long int &  nps 
) [private]

Definition at line 429 of file ReadMCNP5.cpp.

                                                                  {

  // get simulation date and time
  // mcnp   version 5     ld=11242008  probid =  03/23/09 13:38:56
  char line[100];
  file.getline(line, 100);
  date_and_time = line;
  if (debug) std::cout << "date_and_time=| " << date_and_time << std::endl;

  // get simulation title
  // iter Module 4                                                                   
  file.getline(line, 100);
  title = line;
  if (debug) std::cout << "title=| " << title  << std::endl;

  // get number of histories
  // Number of histories used for normalizing tallies =      50000000.00
  file.getline(line, 100);
  std::string a = line;
  std::string::size_type b = a.find("Number of histories used for normalizing tallies =");
  if (std::string::npos!=b) {
    std::istringstream nps_ss( 
      a.substr(b+sizeof("Number of histories used for normalizing tallies ="),100) );
    nps_ss >> nps;
    if (debug) std::cout << "nps=| " << nps << std::endl;
  } else return MB_FAILURE;

  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::read_mesh_planes ( std::fstream &  file,
bool  debug,
std::vector< double >  planes[3],
coordinate_system coord_sys 
) [private]

Definition at line 540 of file ReadMCNP5.cpp.

                                                                          {

  // Tally bin boundaries:
  ErrorCode result;
  char line[10000];
  file.getline(line, 10000);
  std::string a = line;
  if (std::string::npos == a.find("Tally bin boundaries:"))
    return MB_FAILURE;
 
  // decide what coordinate system the tally is using
  // first check for Cylindrical coordinates:
  file.getline(line, 10000);
  a = line;
  std::string::size_type b = a.find("Cylinder origin at");
  if (std::string::npos != b) {
    coord_sys = CYLINDRICAL;
    if (debug) std::cout << "origin, axis, direction=| " << a << std::endl;
    std::istringstream ss(a.substr(b+sizeof("Cylinder origin at"),10000));
    // The meshtal file does not contain sufficient information to transform
    // the particle. Although origin, axs, and vec is needed only origin and
    // axs appear in the meshtal file. Why was vec omitted?.
    // get origin (not used)
    // Cylinder origin at   0.00E+00  0.00E+00  0.00E+00, 
    // axis in  0.000E+00 0.000E+00 1.000E+00 direction
    double origin[3];
    if (debug) std::cout << "origin=| ";
    for (int i=0; i<3; i++) {
      ss >> origin[i];
      if (debug) std::cout << origin[i] << " ";
    }
    if (debug) std::cout << std::endl;
    int length_of_string = 10;
    ss.ignore( length_of_string, ' ');
    ss.ignore( length_of_string, ' ');
    ss.ignore( length_of_string, ' ');
    // get axis (not used)
    double axis[3];
    if (debug) std::cout << "axis=| ";
    for (int i=0; i<3; i++) {
      ss >> axis[i];
      if (debug) std::cout << axis[i] << " ";
    }
    if (debug) std::cout << std::endl;
    file.getline(line, 10000);
    a = line;

    // get r planes
    if (debug) std::cout << "R direction:=| ";
    b = a.find("R direction:");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("R direction"),10000));
      result = get_mesh_plane( ss2, debug, planes[0] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;

    // get z planes
    file.getline(line, 10000);
    a = line;
    if (debug) std::cout << "Z direction:=| ";
    b = a.find("Z direction:");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("Z direction"),10000));
      result = get_mesh_plane( ss2, debug, planes[1] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;

    // get theta planes
    file.getline(line, 10000);
    a = line;
    if (debug) std::cout << "Theta direction:=| ";
    b = a.find("Theta direction (revolutions):");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("Theta direction (revolutions):"),10000));
      result = get_mesh_plane( ss2, debug, planes[2] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;
    
  // Cartesian coordinate system:
  }  else if (std::string::npos != a.find("X direction:") ) {
    coord_sys = CARTESIAN;
    // get x planes
    if (debug) std::cout << "X direction:=| ";
    b = a.find("X direction:");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("X direction"),10000));
      result = get_mesh_plane( ss2, debug, planes[0] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;

    // get y planes
    file.getline(line, 10000);
    a = line;
    if (debug) std::cout << "Y direction:=| ";
    b = a.find("Y direction:");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("Y direction"),10000));
      result = get_mesh_plane( ss2, debug, planes[1] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;

    // get z planes
    file.getline(line, 10000);
    a = line;
    if (debug) std::cout << "Z direction:=| ";
    b = a.find("Z direction:");
    if (std::string::npos != b) {
      std::istringstream ss2(a.substr(b+sizeof("Z direction"),10000));
      result = get_mesh_plane( ss2, debug, planes[2] );
      if(MB_SUCCESS != result) return result;
    } else return MB_FAILURE;

  // Spherical coordinate system not yet implemented:
  } else return MB_FAILURE;
    
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::read_tag_values ( const char *  file_name,
const char *  tag_name,
const FileOptions opts,
std::vector< int > &  tag_values_out,
const SubsetList subset_list = 0 
) [virtual]

Read tag values from a file.

Read the list if all integer tag values from the file for a tag that is a single integer value per entity.

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 59 of file ReadMCNP5.cpp.

{
  return MB_NOT_IMPLEMENTED;
}
ErrorCode moab::ReadMCNP5::read_tally_header ( std::fstream &  file,
bool  debug,
unsigned int &  tally_number,
char  tally_comment[100],
particle tally_particle 
) [private]

Definition at line 480 of file ReadMCNP5.cpp.

                                                                           {

  // get tally number
  // Mesh Tally Number 104
  ErrorCode result;
  char line[100];
  file.getline(line, 100);
  std::string a = line;
  std::string::size_type b = a.find("Mesh Tally Number");
  if (std::string::npos != b) {
    std::istringstream tally_number_ss( a.substr(b+sizeof("Mesh Tally Number"),100) );
    tally_number_ss >> tally_number;
    if (debug) std::cout << "tally_number=| " << tally_number << std::endl;
  } else {
    std::cout << "tally number not found" << std::endl;
    return MB_FAILURE;
  }
 
  // next get the tally comment (optional) and particle type
  // 3mm neutron heating in Be (W/cc)     
  // This is a neutron mesh tally.
  // std::string tally_comment;
  
  // get tally particle
  file.getline(line, 100);
  a = line;
  result = get_tally_particle(a, debug, tally_particle);
  if (MB_FAILURE == result) {
    // If this line does not specify the particle type, then it is a tally comment.
    // Get the comment, then get the particle type from the next line.
    tally_comment = line;
    file.getline(line, 100);
    a = line;
    result = get_tally_particle(a, debug, tally_particle);
    if(MB_SUCCESS != result) return result;
  }
  if (debug) std::cout << "tally_comment=| " << tally_comment << std::endl;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::set_header_tags ( EntityHandle  output_meshset,
char  date_and_time[100],
char  title[100],
unsigned long int  nps,
Tag  data_and_time_tag,
Tag  title_tag,
Tag  nps_tag 
) [private]

Definition at line 463 of file ReadMCNP5.cpp.

                                                                  {
  ErrorCode result;
  result = MBI->tag_set_data( data_and_time_tag, &output_meshset, 1, &date_and_time);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( title_tag, &output_meshset, 1, &title);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( nps_tag, &output_meshset, 1, &nps);
  if(MB_SUCCESS != result) return result;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::set_tally_tags ( EntityHandle  tally_meshset,
unsigned int  tally_number,
char  tally_comment[100],
particle  tally_particle,
coordinate_system  tally_coord_sys,
Tag  tally_number_tag,
Tag  tally_comment_tag,
Tag  tally_particle_tag,
Tag  tally_coord_sys_tag 
) [private]

Definition at line 727 of file ReadMCNP5.cpp.

                                                                             {
  ErrorCode result;
  result = MBI->tag_set_data( tally_number_tag,    &tally_meshset, 1, &tally_number);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( tally_comment_tag,   &tally_meshset, 1, &tally_comment);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( tally_particle_tag,  &tally_meshset, 1, &tally_particle);
  if(MB_SUCCESS != result) return result;
  result = MBI->tag_set_data( tally_coord_sys_tag, &tally_meshset, 1, &tally_coord_sys);
  if(MB_SUCCESS != result) return result;
  return MB_SUCCESS;
}
ErrorCode moab::ReadMCNP5::transform_point_to_cartesian ( double *  in,
double *  out,
coordinate_system  coord_sys 
) [private]

Definition at line 957 of file ReadMCNP5.cpp.

                                                                                  {
  // Transform coordinate system
  switch(coord_sys) {
    case CARTESIAN :
      out[0] = in[0]; // x
      out[1] = in[1]; // y
      out[2] = in[2]; // z
      break;
    // theta is in rotations
    case CYLINDRICAL :
      out[0] = in[0]*cos( 2*PI*in[2] ); // x
      out[1] = in[0]*sin( 2*PI*in[2] ); // y
      out[2] = in[1];                  // z
      break;
    case SPHERICAL :
      return MB_NOT_IMPLEMENTED;
    default :
      return MB_NOT_IMPLEMENTED;
  }
  
  return MB_SUCCESS;
}

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 99 of file ReadMCNP5.hpp.

const Tag* moab::ReadMCNP5::fileIDTag [private]

Definition at line 98 of file ReadMCNP5.hpp.

Definition at line 96 of file ReadMCNP5.hpp.

int moab::ReadMCNP5::nodeId [private]

Definition at line 99 of file ReadMCNP5.hpp.

const double moab::ReadMCNP5::PI = 3.141592653589793 [static, private]

Definition at line 80 of file ReadMCNP5.hpp.

Definition at line 93 of file ReadMCNP5.hpp.


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