MOAB: Mesh Oriented datABase  (version 5.3.0)
ReadMCNP5.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 #include "ReadMCNP5.hpp"
00017 #include "moab/Interface.hpp"
00018 #include "moab/ReadUtilIface.hpp"
00019 #include "Internals.hpp"  // For MB_START_ID
00020 #include "moab/Range.hpp"
00021 #include "moab/FileOptions.hpp"
00022 #include "moab/Util.hpp"
00023 
00024 #include <iostream>
00025 #include <sstream>
00026 #include <fstream>
00027 #include <vector>
00028 #include <cstdlib>
00029 #include <cmath>
00030 #include <cassert>
00031 
00032 namespace moab
00033 {
00034 
00035 // Parameters
00036 const double ReadMCNP5::PI   = 3.141592653589793;
00037 const double ReadMCNP5::C2PI = 0.1591549430918954;
00038 const double ReadMCNP5::CPI  = 0.3183098861837907;
00039 
00040 ReaderIface* ReadMCNP5::factory( Interface* iface )
00041 {
00042     return new ReadMCNP5( iface );
00043 }
00044 
00045 // Constructor
00046 ReadMCNP5::ReadMCNP5( Interface* impl ) : MBI( impl ), fileIDTag( NULL ), nodeId( 0 ), elemId( 0 )
00047 {
00048     assert( NULL != impl );
00049     MBI->query_interface( readMeshIface );
00050     assert( NULL != readMeshIface );
00051 }
00052 
00053 // Destructor
00054 ReadMCNP5::~ReadMCNP5()
00055 {
00056     if( readMeshIface )
00057     {
00058         MBI->release_interface( readMeshIface );
00059         readMeshIface = 0;
00060     }
00061 }
00062 
00063 ErrorCode ReadMCNP5::read_tag_values( const char* /* file_name */, const char* /* tag_name */,
00064                                       const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */,
00065                                       const SubsetList* /* subset_list */ )
00066 {
00067     return MB_NOT_IMPLEMENTED;
00068 }
00069 
00070 // Load the file as called by the Interface function
00071 ErrorCode ReadMCNP5::load_file( const char* filename, const EntityHandle* input_meshset, const FileOptions& options,
00072                                 const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag )
00073 {
00074     // At this time there is no support for reading a subset of the file
00075     if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for meshtal" ); }
00076 
00077     nodeId = elemId = 0;
00078     fileIDTag       = file_id_tag;
00079 
00080     // Average several meshtal files if the AVERAGE_TALLY option is given.
00081     // In this case, the integer value is the number of files to average.
00082     // If averaging, the filename passed to load_file is assumed to be the
00083     // root filename. The files are indexed as "root_filename""index".meshtal.
00084     // Indices start with 1.
00085     int n_files;
00086     bool average = false;
00087     ErrorCode result;
00088     if( MB_SUCCESS == options.get_int_option( "AVERAGE_TALLY", n_files ) )
00089     {
00090         // Read the first file (but do not average -> can't average a single file)
00091         result = load_one_file( filename, input_meshset, options, average );
00092         if( MB_SUCCESS != result ) return result;
00093 
00094         // Get the root filename
00095         std::string root_filename( filename );
00096         int length = root_filename.length();
00097         root_filename.erase( length - sizeof( ".meshtal" ) );
00098 
00099         // Average the first file with the rest of the files
00100         average = true;
00101         for( int i = 2; i <= n_files; i++ )
00102         {
00103             std::stringstream index;
00104             index << i;
00105             std::string subsequent_filename = root_filename + index.str() + ".meshtal";
00106             result = load_one_file( subsequent_filename.c_str(), input_meshset, options, average );
00107             if( MB_SUCCESS != result ) return result;
00108         }
00109 
00110         // If not averaging, read a single file
00111     }
00112     else
00113     {
00114         result = load_one_file( filename, input_meshset, options, average );
00115         if( MB_SUCCESS != result ) return result;
00116     }
00117 
00118     return MB_SUCCESS;
00119 }
00120 
00121 // This actually reads the file. It creates the mesh elements unless
00122 // the file is being averaged with a pre-existing mesh.
00123 ErrorCode ReadMCNP5::load_one_file( const char* fname, const EntityHandle* input_meshset, const FileOptions& options,
00124                                     const bool average )
00125 {
00126     bool debug = false;
00127     if( debug ) std::cout << "begin ReadMCNP5::load_one_file" << std::endl;
00128 
00129     ErrorCode result;
00130     std::fstream file;
00131     file.open( fname, std::fstream::in );
00132     char line[10000];
00133 
00134     // Create tags
00135     Tag date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag, tally_particle_tag,
00136         tally_coord_sys_tag, tally_tag, error_tag;
00137 
00138     result = create_tags( date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag,
00139                           tally_particle_tag, tally_coord_sys_tag, tally_tag, error_tag );
00140     if( MB_SUCCESS != result ) return result;
00141 
00142     // ******************************************************************
00143     // This info exists only at the top of each meshtal file
00144     // ******************************************************************
00145 
00146     // Define characteristics of the entire file
00147     char date_and_time[100] = "";
00148     char title[100]         = "";
00149     // This file's number of particles
00150     unsigned long int nps;
00151     // Sum of this file's and existing file's nps for averaging
00152     unsigned long int new_nps;
00153 
00154     // Read the file header
00155     result = read_file_header( file, debug, date_and_time, title, nps );
00156     if( MB_SUCCESS != result ) return result;
00157 
00158     // Blank line
00159     file.getline( line, 10000 );
00160 
00161     // Everything stored in the file being read will be in the input_meshset.
00162     // If this is being saved in MOAB, set header tags
00163     if( !average && 0 != input_meshset )
00164     {
00165         result = set_header_tags( *input_meshset, date_and_time, title, nps, date_and_time_tag, title_tag, nps_tag );
00166         if( MB_SUCCESS != result ) return result;
00167     }
00168 
00169     // ******************************************************************
00170     // This info is repeated for each tally in the meshtal file.
00171     // ******************************************************************
00172 
00173     // If averaging, nps will hold the sum of particles simulated in both tallies.
00174     while( !file.eof() )
00175     {
00176         // Define characteristics of this tally
00177         unsigned int tally_number;
00178         char tally_comment[100] = "";
00179         particle tally_particle;
00180         coordinate_system tally_coord_sys;
00181         std::vector< double > planes[3];
00182         unsigned int n_chopped_x0_planes;
00183         unsigned int n_chopped_x2_planes;
00184 
00185         // Read tally header
00186         result = read_tally_header( file, debug, tally_number, tally_comment, tally_particle );
00187         if( MB_SUCCESS != result ) return result;
00188 
00189         // Blank line
00190         file.getline( line, 10000 );
00191         std::string l = line;
00192         // if this string is present then skip the following blank line
00193         if( std::string::npos != l.find( "This mesh tally is modified by a dose response function." ) )
00194         {
00195             file.getline( line, 10000 );
00196         }
00197 
00198         // Read mesh planes
00199         result = read_mesh_planes( file, debug, planes, tally_coord_sys );
00200         if( MB_SUCCESS != result ) return result;
00201 
00202         // Get energy boundaries
00203         file.getline( line, 10000 );
00204         std::string a = line;
00205         if( debug ) std::cout << "Energy bin boundaries:=| " << a << std::endl;
00206 
00207         // Blank
00208         file.getline( line, 10000 );
00209 
00210         // Column headers
00211         file.getline( line, 10000 );
00212 
00213         // If using cylindrical mesh, it may be necessary to chop off the last theta element.
00214         // We chop off the last theta plane because the elements will be wrong and skew up
00215         // the tree building code. This is
00216         // because the hex elements are a linear approximation to the cylindrical elements.
00217         // Chopping off the last plane is problem-dependent, and due to MCNP5's mandate
00218         // that the cylindrical mesh must span 360 degrees.
00219         if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_LAST_AZIMUTHAL_PLANE" ) )
00220         {
00221             planes[2].pop_back();
00222             n_chopped_x2_planes = 1;
00223             if( debug ) std::cout << "remove last cylindrical plane option found" << std::endl;
00224         }
00225         else
00226         {
00227             n_chopped_x2_planes = 0;
00228         }
00229 
00230         // If using cylindrical mesh, it may be necessary to chop off the first radial element.
00231         // These elements extend from the axis and often do not cover areas of interest. For
00232         // example, if the mesh was meant to cover r=390-400, the first layer will go from
00233         // 0-390 and serve as incorrect source elements for interpolation.
00234         if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_FIRST_RADIAL_PLANE" ) )
00235         {
00236             std::vector< double >::iterator front = planes[0].begin();
00237             planes[0].erase( front );
00238             n_chopped_x0_planes = 1;
00239             if( debug ) std::cout << "remove first radial plane option found" << std::endl;
00240         }
00241         else
00242         {
00243             n_chopped_x0_planes = 0;
00244         }
00245 
00246         // Read the values and errors of each element from the file.
00247         // Do not read values that are chopped off.
00248         unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
00249         double* values          = new double[n_elements];
00250         double* errors          = new double[n_elements];
00251         result = read_element_values_and_errors( file, debug, planes, n_chopped_x0_planes, n_chopped_x2_planes,
00252                                                  tally_particle, values, errors );
00253         if( MB_SUCCESS != result ) return result;
00254 
00255         // Blank line
00256         file.getline( line, 10000 );
00257 
00258         // ****************************************************************
00259         // This tally has been read. If it is not being averaged, build tags,
00260         // vertices and elements. If it is being averaged, average the data
00261         // with a tally already existing in the MOAB instance.
00262         // ****************************************************************
00263         if( !average )
00264         {
00265             EntityHandle tally_meshset;
00266             result = MBI->create_meshset( MESHSET_SET, tally_meshset );
00267             if( MB_SUCCESS != result ) return result;
00268 
00269             // Set tags on the tally
00270             result = set_tally_tags( tally_meshset, tally_number, tally_comment, tally_particle, tally_coord_sys,
00271                                      tally_number_tag, tally_comment_tag, tally_particle_tag, tally_coord_sys_tag );
00272             if( MB_SUCCESS != result ) return result;
00273 
00274             // The only info needed to build elements is the mesh plane boundaries.
00275             // Build vertices...
00276             EntityHandle start_vert = 0;
00277             result                  = create_vertices( planes, debug, start_vert, tally_coord_sys, tally_meshset );
00278             if( MB_SUCCESS != result ) return result;
00279 
00280             // Build elements and tag them with tally values and errors, then add
00281             // them to the tally_meshset.
00282             result = create_elements( debug, planes, n_chopped_x0_planes, n_chopped_x2_planes, start_vert, values,
00283                                       errors, tally_tag, error_tag, tally_meshset, tally_coord_sys );
00284             if( MB_SUCCESS != result ) return result;
00285 
00286             // Add this tally's meshset to the output meshset
00287             if( debug ) std::cout << "not averaging tally" << std::endl;
00288 
00289             // Average the tally values, then delete stuff that was created
00290         }
00291         else
00292         {
00293             if( debug ) std::cout << "averaging tally" << std::endl;
00294             result = average_with_existing_tally( debug, new_nps, nps, tally_number, tally_number_tag, nps_tag,
00295                                                   tally_tag, error_tag, values, errors, n_elements );
00296             if( MB_SUCCESS != result ) return result;
00297         }
00298 
00299         // Clean up
00300         delete[] values;
00301         delete[] errors;
00302     }
00303 
00304     // If we are averaging, delete the remainder of this file's information.
00305     // Add the new nps to the existing file's nps if we are averaging.
00306     // This is calculated during every tally averaging but only used after the last one.
00307     if( average )
00308     {
00309         Range matching_nps_sets;
00310         result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, matching_nps_sets );
00311         if( MB_SUCCESS != result ) return result;
00312         if( debug ) std::cout << "number of matching nps meshsets=" << matching_nps_sets.size() << std::endl;
00313         assert( 1 == matching_nps_sets.size() );
00314         result = MBI->tag_set_data( nps_tag, matching_nps_sets, &new_nps );
00315         if( MB_SUCCESS != result ) return result;
00316 
00317         // If this file is not being averaged, return the output meshset.
00318     }
00319 
00320     file.close();
00321     return MB_SUCCESS;
00322 }
00323 
00324 // create tags needed for this reader
00325 ErrorCode ReadMCNP5::create_tags( Tag& date_and_time_tag, Tag& title_tag, Tag& nps_tag, Tag& tally_number_tag,
00326                                   Tag& tally_comment_tag, Tag& tally_particle_tag, Tag& tally_coord_sys_tag,
00327                                   Tag& tally_tag, Tag& error_tag )
00328 {
00329     ErrorCode result;
00330     result = MBI->tag_get_handle( "DATE_AND_TIME_TAG", 100, MB_TYPE_OPAQUE, date_and_time_tag,
00331                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00332     if( MB_SUCCESS != result ) return result;
00333     result = MBI->tag_get_handle( "TITLE_TAG", 100, MB_TYPE_OPAQUE, title_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00334     if( MB_SUCCESS != result ) return result;
00335     result = MBI->tag_get_handle( "NPS_TAG", sizeof( unsigned long int ), MB_TYPE_OPAQUE, nps_tag,
00336                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00337     if( MB_SUCCESS != result ) return result;
00338     result =
00339         MBI->tag_get_handle( "TALLY_NUMBER_TAG", 1, MB_TYPE_INTEGER, tally_number_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00340     if( MB_SUCCESS != result ) return result;
00341     result = MBI->tag_get_handle( "TALLY_COMMENT_TAG", 100, MB_TYPE_OPAQUE, tally_comment_tag,
00342                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00343     if( MB_SUCCESS != result ) return result;
00344     result = MBI->tag_get_handle( "TALLY_PARTICLE_TAG", sizeof( particle ), MB_TYPE_OPAQUE, tally_particle_tag,
00345                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00346     if( MB_SUCCESS != result ) return result;
00347     result = MBI->tag_get_handle( "TALLY_COORD_SYS_TAG", sizeof( coordinate_system ), MB_TYPE_OPAQUE,
00348                                   tally_coord_sys_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00349     if( MB_SUCCESS != result ) return result;
00350     result = MBI->tag_get_handle( "TALLY_TAG", 1, MB_TYPE_DOUBLE, tally_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00351     if( MB_SUCCESS != result ) return result;
00352     result = MBI->tag_get_handle( "ERROR_TAG", 1, MB_TYPE_DOUBLE, error_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00353     if( MB_SUCCESS != result ) return result;
00354 
00355     return MB_SUCCESS;
00356 }
00357 
00358 ErrorCode ReadMCNP5::read_file_header( std::fstream& file, bool debug, char date_and_time[100], char title[100],
00359                                        unsigned long int& nps )
00360 {
00361     // Get simulation date and time
00362     // mcnp   version 5     ld=11242008  probid =  03/23/09 13:38:56
00363     char line[100];
00364     file.getline( line, 100 );
00365     date_and_time = line;
00366     if( debug ) std::cout << "date_and_time=| " << date_and_time << std::endl;
00367 
00368     // Get simulation title
00369     // iter Module 4
00370     file.getline( line, 100 );
00371     title = line;
00372     if( debug ) std::cout << "title=| " << title << std::endl;
00373 
00374     // Get number of histories
00375     // Number of histories used for normalizing tallies =      50000000.00
00376     file.getline( line, 100 );
00377     std::string a            = line;
00378     std::string::size_type b = a.find( "Number of histories used for normalizing tallies =" );
00379     if( std::string::npos != b )
00380     {
00381         std::istringstream nps_ss(
00382             a.substr( b + sizeof( "Number of histories used for normalizing tallies =" ), 100 ) );
00383         nps_ss >> nps;
00384         if( debug ) std::cout << "nps=| " << nps << std::endl;
00385     }
00386     else
00387         return MB_FAILURE;
00388 
00389     return MB_SUCCESS;
00390 }
00391 
00392 ErrorCode ReadMCNP5::set_header_tags( EntityHandle output_meshset, char date_and_time[100], char title[100],
00393                                       unsigned long int nps, Tag data_and_time_tag, Tag title_tag, Tag nps_tag )
00394 {
00395     ErrorCode result;
00396     result = MBI->tag_set_data( data_and_time_tag, &output_meshset, 1, &date_and_time );
00397     if( MB_SUCCESS != result ) return result;
00398     result = MBI->tag_set_data( title_tag, &output_meshset, 1, &title );
00399     if( MB_SUCCESS != result ) return result;
00400     result = MBI->tag_set_data( nps_tag, &output_meshset, 1, &nps );
00401     if( MB_SUCCESS != result ) return result;
00402 
00403     return MB_SUCCESS;
00404 }
00405 
00406 ErrorCode ReadMCNP5::read_tally_header( std::fstream& file, bool debug, unsigned int& tally_number,
00407                                         char tally_comment[100], particle& tally_particle )
00408 {
00409     // Get tally number
00410     // Mesh Tally Number 104
00411     ErrorCode result;
00412     char line[100];
00413     file.getline( line, 100 );
00414     std::string a            = line;
00415     std::string::size_type b = a.find( "Mesh Tally Number" );
00416     if( std::string::npos != b )
00417     {
00418         std::istringstream tally_number_ss( a.substr( b + sizeof( "Mesh Tally Number" ), 100 ) );
00419         tally_number_ss >> tally_number;
00420         if( debug ) std::cout << "tally_number=| " << tally_number << std::endl;
00421     }
00422     else
00423     {
00424         std::cout << "tally number not found" << std::endl;
00425         return MB_FAILURE;
00426     }
00427 
00428     // Next get the tally comment (optional) and particle type
00429     // 3mm neutron heating in Be (W/cc)
00430     // This is a neutron mesh tally.
00431     // std::string tally_comment;
00432 
00433     // Get tally particle
00434     file.getline( line, 100 );
00435     a      = line;
00436     result = get_tally_particle( a, debug, tally_particle );
00437     if( MB_FAILURE == result )
00438     {
00439         // If this line does not specify the particle type, then it is a tally comment.
00440         // Get the comment, then get the particle type from the next line.
00441         tally_comment = line;
00442         file.getline( line, 100 );
00443         a      = line;
00444         result = get_tally_particle( a, debug, tally_particle );
00445         if( MB_SUCCESS != result ) return result;
00446     }
00447     if( debug ) std::cout << "tally_comment=| " << tally_comment << std::endl;
00448 
00449     return MB_SUCCESS;
00450 }
00451 
00452 ErrorCode ReadMCNP5::get_tally_particle( std::string a, bool debug, particle& tally_particle )
00453 {
00454     if( std::string::npos != a.find( "This is a neutron mesh tally." ) ) { tally_particle = NEUTRON; }
00455     else if( std::string::npos != a.find( "This is a photon mesh tally." ) )
00456     {
00457         tally_particle = PHOTON;
00458     }
00459     else if( std::string::npos != a.find( "This is an electron mesh tally." ) )
00460     {
00461         tally_particle = ELECTRON;
00462     }
00463     else
00464         return MB_FAILURE;
00465 
00466     if( debug ) std::cout << "tally_particle=| " << tally_particle << std::endl;
00467 
00468     return MB_SUCCESS;
00469 }
00470 
00471 ErrorCode ReadMCNP5::read_mesh_planes( std::fstream& file, bool debug, std::vector< double > planes[3],
00472                                        coordinate_system& coord_sys )
00473 {
00474     // Tally bin boundaries:
00475     ErrorCode result;
00476     char line[10000];
00477     file.getline( line, 10000 );
00478     std::string a = line;
00479     if( std::string::npos == a.find( "Tally bin boundaries:" ) ) return MB_FAILURE;
00480 
00481     // Decide what coordinate system the tally is using
00482     // first check for Cylindrical coordinates:
00483     file.getline( line, 10000 );
00484     a                        = line;
00485     std::string::size_type b = a.find( "Cylinder origin at" );
00486     if( std::string::npos != b )
00487     {
00488         coord_sys = CYLINDRICAL;
00489         if( debug ) std::cout << "origin, axis, direction=| " << a << std::endl;
00490         std::istringstream ss( a.substr( b + sizeof( "Cylinder origin at" ), 10000 ) );
00491         // The meshtal file does not contain sufficient information to transform
00492         // the particle. Although origin, axs, and vec is needed only origin and
00493         // axs appear in the meshtal file. Why was vec omitted?.
00494         // get origin (not used)
00495         // Cylinder origin at   0.00E+00  0.00E+00  0.00E+00,
00496         // axis in  0.000E+00 0.000E+00 1.000E+00 direction
00497         double origin[3];
00498         if( debug ) std::cout << "origin=| ";
00499         for( int i = 0; i < 3; i++ )
00500         {
00501             ss >> origin[i];
00502             if( debug ) std::cout << origin[i] << " ";
00503         }
00504         if( debug ) std::cout << std::endl;
00505         int length_of_string = 10;
00506         ss.ignore( length_of_string, ' ' );
00507         ss.ignore( length_of_string, ' ' );
00508         ss.ignore( length_of_string, ' ' );
00509         // Get axis (not used)
00510         double axis[3];
00511         if( debug ) std::cout << "axis=| ";
00512         for( int i = 0; i < 3; i++ )
00513         {
00514             ss >> axis[i];
00515             if( debug ) std::cout << axis[i] << " ";
00516         }
00517         if( debug ) std::cout << std::endl;
00518         file.getline( line, 10000 );
00519         a = line;
00520 
00521         // Get r planes
00522         if( debug ) std::cout << "R direction:=| ";
00523         b = a.find( "R direction:" );
00524         if( std::string::npos != b )
00525         {
00526             std::istringstream ss2( a.substr( b + sizeof( "R direction" ), 10000 ) );
00527             result = get_mesh_plane( ss2, debug, planes[0] );
00528             if( MB_SUCCESS != result ) return result;
00529         }
00530         else
00531             return MB_FAILURE;
00532 
00533         // Get z planes
00534         file.getline( line, 10000 );
00535         a = line;
00536         if( debug ) std::cout << "Z direction:=| ";
00537         b = a.find( "Z direction:" );
00538         if( std::string::npos != b )
00539         {
00540             std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
00541             result = get_mesh_plane( ss2, debug, planes[1] );
00542             if( MB_SUCCESS != result ) return result;
00543         }
00544         else
00545             return MB_FAILURE;
00546 
00547         // Get theta planes
00548         file.getline( line, 10000 );
00549         a = line;
00550         if( debug ) std::cout << "Theta direction:=| ";
00551         b = a.find( "Theta direction (revolutions):" );
00552         if( std::string::npos != b )
00553         {
00554             std::istringstream ss2( a.substr( b + sizeof( "Theta direction (revolutions):" ), 10000 ) );
00555             result = get_mesh_plane( ss2, debug, planes[2] );
00556             if( MB_SUCCESS != result ) return result;
00557         }
00558         else
00559             return MB_FAILURE;
00560 
00561         // Cartesian coordinate system:
00562     }
00563     else if( std::string::npos != a.find( "X direction:" ) )
00564     {
00565         coord_sys = CARTESIAN;
00566         // Get x planes
00567         if( debug ) std::cout << "X direction:=| ";
00568         b = a.find( "X direction:" );
00569         if( std::string::npos != b )
00570         {
00571             std::istringstream ss2( a.substr( b + sizeof( "X direction" ), 10000 ) );
00572             result = get_mesh_plane( ss2, debug, planes[0] );
00573             if( MB_SUCCESS != result ) return result;
00574         }
00575         else
00576             return MB_FAILURE;
00577 
00578         // Get y planes
00579         file.getline( line, 10000 );
00580         a = line;
00581         if( debug ) std::cout << "Y direction:=| ";
00582         b = a.find( "Y direction:" );
00583         if( std::string::npos != b )
00584         {
00585             std::istringstream ss2( a.substr( b + sizeof( "Y direction" ), 10000 ) );
00586             result = get_mesh_plane( ss2, debug, planes[1] );
00587             if( MB_SUCCESS != result ) return result;
00588         }
00589         else
00590             return MB_FAILURE;
00591 
00592         // Get z planes
00593         file.getline( line, 10000 );
00594         a = line;
00595         if( debug ) std::cout << "Z direction:=| ";
00596         b = a.find( "Z direction:" );
00597         if( std::string::npos != b )
00598         {
00599             std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
00600             result = get_mesh_plane( ss2, debug, planes[2] );
00601             if( MB_SUCCESS != result ) return result;
00602         }
00603         else
00604             return MB_FAILURE;
00605 
00606         // Spherical coordinate system not yet implemented:
00607     }
00608     else
00609         return MB_FAILURE;
00610 
00611     return MB_SUCCESS;
00612 }
00613 
00614 // Given a stringstream, return a vector of values in the string.
00615 ErrorCode ReadMCNP5::get_mesh_plane( std::istringstream& ss, bool debug, std::vector< double >& plane )
00616 {
00617     double value;
00618     plane.clear();
00619     while( !ss.eof() )
00620     {
00621         ss >> value;
00622         plane.push_back( value );
00623         if( debug ) std::cout << value << " ";
00624     }
00625     if( debug ) std::cout << std::endl;
00626 
00627     return MB_SUCCESS;
00628 }
00629 
00630 ErrorCode ReadMCNP5::read_element_values_and_errors( std::fstream& file, bool /* debug */,
00631                                                      std::vector< double > planes[3], unsigned int n_chopped_x0_planes,
00632                                                      unsigned int n_chopped_x2_planes, particle tally_particle,
00633                                                      double values[], double errors[] )
00634 {
00635     unsigned int index = 0;
00636     // Need to read every line in the file, even if we chop off some elements
00637     for( unsigned int i = 0; i < planes[0].size() - 1 + n_chopped_x0_planes; i++ )
00638     {
00639         for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
00640         {
00641             for( unsigned int k = 0; k < planes[2].size() - 1 + n_chopped_x2_planes; k++ )
00642             {
00643                 char line[100];
00644                 file.getline( line, 100 );
00645                 // If this element has been chopped off, skip it
00646                 if( i < n_chopped_x0_planes ) continue;
00647                 if( k >= planes[2].size() - 1 && k < planes[2].size() - 1 + n_chopped_x2_planes ) continue;
00648                 std::string a = line;
00649                 std::stringstream ss( a );
00650                 double centroid[3];
00651                 double energy;
00652                 // For some reason, photon tallies print the energy in the first column
00653                 if( PHOTON == tally_particle ) ss >> energy;
00654                 // The centroid is not used in this reader
00655                 ss >> centroid[0];
00656                 ss >> centroid[1];
00657                 ss >> centroid[2];
00658                 // Only the tally values and errors are used
00659                 ss >> values[index];
00660                 ss >> errors[index];
00661 
00662                 // Make sure that input data is good
00663                 if( !Util::is_finite( errors[index] ) )
00664                 {
00665                     std::cerr << "found nan error while reading file" << std::endl;
00666                     errors[index] = 1.0;
00667                 }
00668                 if( !Util::is_finite( values[index] ) )
00669                 {
00670                     std::cerr << "found nan value while reading file" << std::endl;
00671                     values[index] = 0.0;
00672                 }
00673 
00674                 index++;
00675             }
00676         }
00677     }
00678 
00679     return MB_SUCCESS;
00680 }
00681 
00682 ErrorCode ReadMCNP5::set_tally_tags( EntityHandle tally_meshset, unsigned int tally_number, char tally_comment[100],
00683                                      particle tally_particle, coordinate_system tally_coord_sys, Tag tally_number_tag,
00684                                      Tag tally_comment_tag, Tag tally_particle_tag, Tag tally_coord_sys_tag )
00685 {
00686     ErrorCode result;
00687     result = MBI->tag_set_data( tally_number_tag, &tally_meshset, 1, &tally_number );
00688     if( MB_SUCCESS != result ) return result;
00689     result = MBI->tag_set_data( tally_comment_tag, &tally_meshset, 1, &tally_comment );
00690     if( MB_SUCCESS != result ) return result;
00691     result = MBI->tag_set_data( tally_particle_tag, &tally_meshset, 1, &tally_particle );
00692     if( MB_SUCCESS != result ) return result;
00693     result = MBI->tag_set_data( tally_coord_sys_tag, &tally_meshset, 1, &tally_coord_sys );
00694     if( MB_SUCCESS != result ) return result;
00695 
00696     return MB_SUCCESS;
00697 }
00698 
00699 ErrorCode ReadMCNP5::create_vertices( std::vector< double > planes[3], bool debug, EntityHandle& start_vert,
00700                                       coordinate_system coord_sys, EntityHandle tally_meshset )
00701 {
00702     // The only info needed to build elements is the mesh plane boundaries.
00703     ErrorCode result;
00704     int n_verts = planes[0].size() * planes[1].size() * planes[2].size();
00705     if( debug ) std::cout << "n_verts=" << n_verts << std::endl;
00706     std::vector< double* > coord_arrays( 3 );
00707     result = readMeshIface->get_node_coords( 3, n_verts, MB_START_ID, start_vert, coord_arrays );
00708     if( MB_SUCCESS != result ) return result;
00709     assert( 0 != start_vert );  // Check for NULL
00710 
00711     for( unsigned int k = 0; k < planes[2].size(); k++ )
00712     {
00713         for( unsigned int j = 0; j < planes[1].size(); j++ )
00714         {
00715             for( unsigned int i = 0; i < planes[0].size(); i++ )
00716             {
00717                 unsigned int idx = ( k * planes[0].size() * planes[1].size() + j * planes[0].size() + i );
00718                 double in[3], out[3];
00719 
00720                 in[0]  = planes[0][i];
00721                 in[1]  = planes[1][j];
00722                 in[2]  = planes[2][k];
00723                 result = transform_point_to_cartesian( in, out, coord_sys );
00724                 if( MB_SUCCESS != result ) return result;
00725 
00726                 // Cppcheck warning (false positive): variable coord_arrays is assigned a value that
00727                 // is never used
00728                 coord_arrays[0][idx] = out[0];
00729                 coord_arrays[1][idx] = out[1];
00730                 coord_arrays[2][idx] = out[2];
00731             }
00732         }
00733     }
00734     Range vert_range( start_vert, start_vert + n_verts - 1 );
00735     result = MBI->add_entities( tally_meshset, vert_range );
00736     if( MB_SUCCESS != result ) return result;
00737 
00738     if( fileIDTag )
00739     {
00740         result = readMeshIface->assign_ids( *fileIDTag, vert_range, nodeId );
00741         if( MB_SUCCESS != result ) return result;
00742         nodeId += vert_range.size();
00743     }
00744 
00745     return MB_SUCCESS;
00746 }
00747 
00748 ErrorCode ReadMCNP5::create_elements( bool debug, std::vector< double > planes[3], unsigned int /*n_chopped_x0_planes*/,
00749                                       unsigned int /*n_chopped_x2_planes*/, EntityHandle start_vert, double* values,
00750                                       double* errors, Tag tally_tag, Tag error_tag, EntityHandle tally_meshset,
00751                                       coordinate_system tally_coord_sys )
00752 {
00753     ErrorCode result;
00754     unsigned int index;
00755     EntityHandle start_element = 0;
00756     unsigned int n_elements    = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
00757     EntityHandle* connect;
00758     result = readMeshIface->get_element_connect( n_elements, 8, MBHEX, MB_START_ID, start_element, connect );
00759     if( MB_SUCCESS != result ) return result;
00760     assert( 0 != start_element );  // Check for NULL
00761 
00762     unsigned int counter = 0;
00763     for( unsigned int i = 0; i < planes[0].size() - 1; i++ )
00764     {
00765         for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
00766         {
00767             for( unsigned int k = 0; k < planes[2].size() - 1; k++ )
00768             {
00769                 index = start_vert + i + j * planes[0].size() + k * planes[0].size() * planes[1].size();
00770                 // For rectangular mesh, the file prints: x y z
00771                 // z changes the fastest and x changes the slowest.
00772                 // This means that the connectivity ordering is not consistent between
00773                 // rectangular and cylindrical mesh.
00774                 if( CARTESIAN == tally_coord_sys )
00775                 {
00776                     connect[0] = index;
00777                     connect[1] = index + 1;
00778                     connect[2] = index + 1 + planes[0].size();
00779                     connect[3] = index + planes[0].size();
00780                     connect[4] = index + planes[0].size() * planes[1].size();
00781                     connect[5] = index + 1 + planes[0].size() * planes[1].size();
00782                     connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
00783                     connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
00784                     // For cylindrical mesh, the file prints: r z theta
00785                     // Theta changes the fastest and r changes the slowest.
00786                 }
00787                 else if( CYLINDRICAL == tally_coord_sys )
00788                 {
00789                     connect[0] = index;
00790                     connect[1] = index + 1;
00791                     connect[2] = index + 1 + planes[0].size() * planes[1].size();
00792                     connect[3] = index + planes[0].size() * planes[1].size();
00793                     connect[4] = index + planes[0].size();
00794                     connect[5] = index + 1 + planes[0].size();
00795                     connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
00796                     connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
00797                 }
00798                 else
00799                     return MB_NOT_IMPLEMENTED;
00800 
00801                 connect += 8;
00802                 counter++;
00803             }
00804         }
00805     }
00806     if( counter != n_elements ) std::cout << "counter=" << counter << " n_elements=" << n_elements << std::endl;
00807 
00808     Range element_range( start_element, start_element + n_elements - 1 );
00809     result = MBI->tag_set_data( tally_tag, element_range, values );
00810     if( MB_SUCCESS != result ) return result;
00811     result = MBI->tag_set_data( error_tag, element_range, errors );
00812     if( MB_SUCCESS != result ) return result;
00813 
00814     // Add the elements to the tally set
00815     result = MBI->add_entities( tally_meshset, element_range );
00816     if( MB_SUCCESS != result ) return result;
00817     if( debug ) std::cout << "Read " << n_elements << " elements from tally." << std::endl;
00818 
00819     if( fileIDTag )
00820     {
00821         result = readMeshIface->assign_ids( *fileIDTag, element_range, elemId );
00822         if( MB_SUCCESS != result ) return result;
00823         elemId += element_range.size();
00824     }
00825 
00826     return MB_SUCCESS;
00827 }
00828 
00829 // Average a tally that was recently read in with one that already exists in
00830 // the interface. Only the existing values will be updated.
00831 ErrorCode ReadMCNP5::average_with_existing_tally( bool debug, unsigned long int& new_nps, unsigned long int nps1,
00832                                                   unsigned int tally_number, Tag tally_number_tag, Tag nps_tag,
00833                                                   Tag tally_tag, Tag error_tag, double* values1, double* errors1,
00834                                                   unsigned int n_elements )
00835 {
00836     // Get the tally number
00837     ErrorCode result;
00838 
00839     // Match the tally number with one from the existing meshtal file
00840     Range matching_tally_number_sets;
00841     const void* const tally_number_val[] = { &tally_number };
00842     result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &tally_number_tag, tally_number_val, 1,
00843                                                 matching_tally_number_sets );
00844     if( MB_SUCCESS != result ) return result;
00845     if( debug ) std::cout << "number of matching meshsets=" << matching_tally_number_sets.size() << std::endl;
00846     assert( 1 == matching_tally_number_sets.size() );
00847 
00848     // Identify which of the meshsets is existing
00849     EntityHandle existing_meshset;
00850     existing_meshset = matching_tally_number_sets.front();
00851 
00852     // Get the existing elements from the set
00853     Range existing_elements;
00854     result = MBI->get_entities_by_type( existing_meshset, MBHEX, existing_elements );
00855     if( MB_SUCCESS != result ) return result;
00856     assert( existing_elements.size() == n_elements );
00857 
00858     // Get the nps of the existing and new tally
00859     unsigned long int nps0;
00860     Range sets_with_this_tag;
00861     result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, sets_with_this_tag );
00862     if( MB_SUCCESS != result ) return result;
00863     if( debug ) std::cout << "number of nps sets=" << sets_with_this_tag.size() << std::endl;
00864     assert( 1 == sets_with_this_tag.size() );
00865     result = MBI->tag_get_data( nps_tag, &sets_with_this_tag.front(), 1, &nps0 );
00866     if( MB_SUCCESS != result ) return result;
00867     if( debug ) std::cout << "nps0=" << nps0 << " nps1=" << nps1 << std::endl;
00868     new_nps = nps0 + nps1;
00869 
00870     // Get tally values from the existing elements
00871     double* values0 = new double[existing_elements.size()];
00872     double* errors0 = new double[existing_elements.size()];
00873     result          = MBI->tag_get_data( tally_tag, existing_elements, values0 );
00874     if( MB_SUCCESS != result )
00875     {
00876         delete[] values0;
00877         delete[] errors0;
00878         return result;
00879     }
00880     result = MBI->tag_get_data( error_tag, existing_elements, errors0 );
00881     if( MB_SUCCESS != result )
00882     {
00883         delete[] values0;
00884         delete[] errors0;
00885         return result;
00886     }
00887 
00888     // Average the values and errors
00889     result = average_tally_values( nps0, nps1, values0, values1, errors0, errors1, n_elements );
00890     if( MB_SUCCESS != result )
00891     {
00892         delete[] values0;
00893         delete[] errors0;
00894         return result;
00895     }
00896 
00897     // Set the averaged information back onto the existing elements
00898     result = MBI->tag_set_data( tally_tag, existing_elements, values0 );
00899     if( MB_SUCCESS != result )
00900     {
00901         delete[] values0;
00902         delete[] errors0;
00903         return result;
00904     }
00905     result = MBI->tag_set_data( error_tag, existing_elements, errors0 );
00906     if( MB_SUCCESS != result )
00907     {
00908         delete[] values0;
00909         delete[] errors0;
00910         return result;
00911     }
00912 
00913     // Cleanup
00914     delete[] values0;
00915     delete[] errors0;
00916 
00917     return MB_SUCCESS;
00918 }
00919 
00920 ErrorCode ReadMCNP5::transform_point_to_cartesian( double* in, double* out, coordinate_system coord_sys )
00921 {
00922     // Transform coordinate system
00923     switch( coord_sys )
00924     {
00925         case CARTESIAN:
00926             out[0] = in[0];  // x
00927             out[1] = in[1];  // y
00928             out[2] = in[2];  // z
00929             break;
00930         // theta is in rotations
00931         case CYLINDRICAL:
00932             out[0] = in[0] * cos( 2 * PI * in[2] );  // x
00933             out[1] = in[0] * sin( 2 * PI * in[2] );  // y
00934             out[2] = in[1];                          // z
00935             break;
00936         case SPHERICAL:
00937             return MB_NOT_IMPLEMENTED;
00938         default:
00939             return MB_NOT_IMPLEMENTED;
00940     }
00941 
00942     return MB_SUCCESS;
00943 }
00944 
00945 // Average two tally values and their error. Return average values in the
00946 // place of first tally values.
00947 ErrorCode ReadMCNP5::average_tally_values( const unsigned long int nps0, const unsigned long int nps1, double* values0,
00948                                            const double* values1, double* errors0, const double* errors1,
00949                                            const unsigned long int n_values )
00950 {
00951     for( unsigned long int i = 0; i < n_values; i++ )
00952     {
00953         // std::cout << " values0=" << values0[i] << " values1=" << values1[i]
00954         //          << " errors0=" << errors0[i] << " errors1=" << errors1[i]
00955         //          << " nps0=" << nps0 << " nps1=" << nps1 << std::endl;
00956         errors0[i] = sqrt( pow( values0[i] * errors0[i] * nps0, 2 ) + pow( values1[i] * errors1[i] * nps1, 2 ) ) /
00957                      ( values0[i] * nps0 + values1[i] * nps1 );
00958 
00959         // It is possible to introduce nans if the values are zero.
00960         if( !Util::is_finite( errors0[i] ) ) errors0[i] = 1.0;
00961 
00962         values0[i] = ( values0[i] * nps0 + values1[i] * nps1 ) / ( nps0 + nps1 );
00963 
00964         // std::cout << " values0=" << values0[i] << " errors0=" << errors0[i] << std::endl;
00965     }
00966     // REMEMBER TO UPDATE NPS0 = NPS0 + NPS1 after this
00967 
00968     return MB_SUCCESS;
00969 }
00970 
00971 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines