MOAB: Mesh Oriented datABase  (version 5.2.1)
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         { file.getline( line, 10000 ); }
00195 
00196         // Read mesh planes
00197         result = read_mesh_planes( file, debug, planes, tally_coord_sys );
00198         if( MB_SUCCESS != result ) return result;
00199 
00200         // Get energy boundaries
00201         file.getline( line, 10000 );
00202         std::string a = line;
00203         if( debug ) std::cout << "Energy bin boundaries:=| " << a << std::endl;
00204 
00205         // Blank
00206         file.getline( line, 10000 );
00207 
00208         // Column headers
00209         file.getline( line, 10000 );
00210 
00211         // If using cylindrical mesh, it may be necessary to chop off the last theta element.
00212         // We chop off the last theta plane because the elements will be wrong and skew up
00213         // the tree building code. This is
00214         // because the hex elements are a linear approximation to the cylindrical elements.
00215         // Chopping off the last plane is problem-dependent, and due to MCNP5's mandate
00216         // that the cylindrical mesh must span 360 degrees.
00217         if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_LAST_AZIMUTHAL_PLANE" ) )
00218         {
00219             planes[2].pop_back();
00220             n_chopped_x2_planes = 1;
00221             if( debug ) std::cout << "remove last cylindrical plane option found" << std::endl;
00222         }
00223         else
00224         {
00225             n_chopped_x2_planes = 0;
00226         }
00227 
00228         // If using cylindrical mesh, it may be necessary to chop off the first radial element.
00229         // These elements extend from the axis and often do not cover areas of interest. For
00230         // example, if the mesh was meant to cover r=390-400, the first layer will go from
00231         // 0-390 and serve as incorrect source elements for interpolation.
00232         if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_FIRST_RADIAL_PLANE" ) )
00233         {
00234             std::vector< double >::iterator front = planes[0].begin();
00235             planes[0].erase( front );
00236             n_chopped_x0_planes = 1;
00237             if( debug ) std::cout << "remove first radial plane option found" << std::endl;
00238         }
00239         else
00240         {
00241             n_chopped_x0_planes = 0;
00242         }
00243 
00244         // Read the values and errors of each element from the file.
00245         // Do not read values that are chopped off.
00246         unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
00247         double* values          = new double[n_elements];
00248         double* errors          = new double[n_elements];
00249         result = read_element_values_and_errors( file, debug, planes, n_chopped_x0_planes, n_chopped_x2_planes,
00250                                                  tally_particle, values, errors );
00251         if( MB_SUCCESS != result ) return result;
00252 
00253         // Blank line
00254         file.getline( line, 10000 );
00255 
00256         // ****************************************************************
00257         // This tally has been read. If it is not being averaged, build tags,
00258         // vertices and elements. If it is being averaged, average the data
00259         // with a tally already existing in the MOAB instance.
00260         // ****************************************************************
00261         if( !average )
00262         {
00263             EntityHandle tally_meshset;
00264             result = MBI->create_meshset( MESHSET_SET, tally_meshset );
00265             if( MB_SUCCESS != result ) return result;
00266 
00267             // Set tags on the tally
00268             result = set_tally_tags( tally_meshset, tally_number, tally_comment, tally_particle, tally_coord_sys,
00269                                      tally_number_tag, tally_comment_tag, tally_particle_tag, tally_coord_sys_tag );
00270             if( MB_SUCCESS != result ) return result;
00271 
00272             // The only info needed to build elements is the mesh plane boundaries.
00273             // Build vertices...
00274             EntityHandle start_vert = 0;
00275             result                  = create_vertices( planes, debug, start_vert, tally_coord_sys, tally_meshset );
00276             if( MB_SUCCESS != result ) return result;
00277 
00278             // Build elements and tag them with tally values and errors, then add
00279             // them to the tally_meshset.
00280             result = create_elements( debug, planes, n_chopped_x0_planes, n_chopped_x2_planes, start_vert, values,
00281                                       errors, tally_tag, error_tag, tally_meshset, tally_coord_sys );
00282             if( MB_SUCCESS != result ) return result;
00283 
00284             // Add this tally's meshset to the output meshset
00285             if( debug ) std::cout << "not averaging tally" << std::endl;
00286 
00287             // Average the tally values, then delete stuff that was created
00288         }
00289         else
00290         {
00291             if( debug ) std::cout << "averaging tally" << std::endl;
00292             result = average_with_existing_tally( debug, new_nps, nps, tally_number, tally_number_tag, nps_tag,
00293                                                   tally_tag, error_tag, values, errors, n_elements );
00294             if( MB_SUCCESS != result ) return result;
00295         }
00296 
00297         // Clean up
00298         delete[] values;
00299         delete[] errors;
00300     }
00301 
00302     // If we are averaging, delete the remainder of this file's information.
00303     // Add the new nps to the existing file's nps if we are averaging.
00304     // This is calculated during every tally averaging but only used after the last one.
00305     if( average )
00306     {
00307         Range matching_nps_sets;
00308         result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, matching_nps_sets );
00309         if( MB_SUCCESS != result ) return result;
00310         if( debug ) std::cout << "number of matching nps meshsets=" << matching_nps_sets.size() << std::endl;
00311         assert( 1 == matching_nps_sets.size() );
00312         result = MBI->tag_set_data( nps_tag, matching_nps_sets, &new_nps );
00313         if( MB_SUCCESS != result ) return result;
00314 
00315         // If this file is not being averaged, return the output meshset.
00316     }
00317 
00318     file.close();
00319     return MB_SUCCESS;
00320 }
00321 
00322 // create tags needed for this reader
00323 ErrorCode ReadMCNP5::create_tags( Tag& date_and_time_tag, Tag& title_tag, Tag& nps_tag, Tag& tally_number_tag,
00324                                   Tag& tally_comment_tag, Tag& tally_particle_tag, Tag& tally_coord_sys_tag,
00325                                   Tag& tally_tag, Tag& error_tag )
00326 {
00327     ErrorCode result;
00328     result = MBI->tag_get_handle( "DATE_AND_TIME_TAG", 100, MB_TYPE_OPAQUE, date_and_time_tag,
00329                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00330     if( MB_SUCCESS != result ) return result;
00331     result = MBI->tag_get_handle( "TITLE_TAG", 100, MB_TYPE_OPAQUE, title_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00332     if( MB_SUCCESS != result ) return result;
00333     result = MBI->tag_get_handle( "NPS_TAG", sizeof( unsigned long int ), MB_TYPE_OPAQUE, nps_tag,
00334                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00335     if( MB_SUCCESS != result ) return result;
00336     result =
00337         MBI->tag_get_handle( "TALLY_NUMBER_TAG", 1, MB_TYPE_INTEGER, tally_number_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00338     if( MB_SUCCESS != result ) return result;
00339     result = MBI->tag_get_handle( "TALLY_COMMENT_TAG", 100, MB_TYPE_OPAQUE, tally_comment_tag,
00340                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00341     if( MB_SUCCESS != result ) return result;
00342     result = MBI->tag_get_handle( "TALLY_PARTICLE_TAG", sizeof( particle ), MB_TYPE_OPAQUE, tally_particle_tag,
00343                                   MB_TAG_SPARSE | MB_TAG_CREAT );
00344     if( MB_SUCCESS != result ) return result;
00345     result = MBI->tag_get_handle( "TALLY_COORD_SYS_TAG", sizeof( coordinate_system ), MB_TYPE_OPAQUE,
00346                                   tally_coord_sys_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00347     if( MB_SUCCESS != result ) return result;
00348     result = MBI->tag_get_handle( "TALLY_TAG", 1, MB_TYPE_DOUBLE, tally_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00349     if( MB_SUCCESS != result ) return result;
00350     result = MBI->tag_get_handle( "ERROR_TAG", 1, MB_TYPE_DOUBLE, error_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00351     if( MB_SUCCESS != result ) return result;
00352 
00353     return MB_SUCCESS;
00354 }
00355 
00356 ErrorCode ReadMCNP5::read_file_header( std::fstream& file, bool debug, char date_and_time[100], char title[100],
00357                                        unsigned long int& nps )
00358 {
00359     // Get simulation date and time
00360     // mcnp   version 5     ld=11242008  probid =  03/23/09 13:38:56
00361     char line[100];
00362     file.getline( line, 100 );
00363     date_and_time = line;
00364     if( debug ) std::cout << "date_and_time=| " << date_and_time << std::endl;
00365 
00366     // Get simulation title
00367     // iter Module 4
00368     file.getline( line, 100 );
00369     title = line;
00370     if( debug ) std::cout << "title=| " << title << std::endl;
00371 
00372     // Get number of histories
00373     // Number of histories used for normalizing tallies =      50000000.00
00374     file.getline( line, 100 );
00375     std::string a            = line;
00376     std::string::size_type b = a.find( "Number of histories used for normalizing tallies =" );
00377     if( std::string::npos != b )
00378     {
00379         std::istringstream nps_ss(
00380             a.substr( b + sizeof( "Number of histories used for normalizing tallies =" ), 100 ) );
00381         nps_ss >> nps;
00382         if( debug ) std::cout << "nps=| " << nps << std::endl;
00383     }
00384     else
00385         return MB_FAILURE;
00386 
00387     return MB_SUCCESS;
00388 }
00389 
00390 ErrorCode ReadMCNP5::set_header_tags( EntityHandle output_meshset, char date_and_time[100], char title[100],
00391                                       unsigned long int nps, Tag data_and_time_tag, Tag title_tag, Tag nps_tag )
00392 {
00393     ErrorCode result;
00394     result = MBI->tag_set_data( data_and_time_tag, &output_meshset, 1, &date_and_time );
00395     if( MB_SUCCESS != result ) return result;
00396     result = MBI->tag_set_data( title_tag, &output_meshset, 1, &title );
00397     if( MB_SUCCESS != result ) return result;
00398     result = MBI->tag_set_data( nps_tag, &output_meshset, 1, &nps );
00399     if( MB_SUCCESS != result ) return result;
00400 
00401     return MB_SUCCESS;
00402 }
00403 
00404 ErrorCode ReadMCNP5::read_tally_header( std::fstream& file, bool debug, unsigned int& tally_number,
00405                                         char tally_comment[100], particle& tally_particle )
00406 {
00407     // Get tally number
00408     // Mesh Tally Number 104
00409     ErrorCode result;
00410     char line[100];
00411     file.getline( line, 100 );
00412     std::string a            = line;
00413     std::string::size_type b = a.find( "Mesh Tally Number" );
00414     if( std::string::npos != b )
00415     {
00416         std::istringstream tally_number_ss( a.substr( b + sizeof( "Mesh Tally Number" ), 100 ) );
00417         tally_number_ss >> tally_number;
00418         if( debug ) std::cout << "tally_number=| " << tally_number << std::endl;
00419     }
00420     else
00421     {
00422         std::cout << "tally number not found" << std::endl;
00423         return MB_FAILURE;
00424     }
00425 
00426     // Next get the tally comment (optional) and particle type
00427     // 3mm neutron heating in Be (W/cc)
00428     // This is a neutron mesh tally.
00429     // std::string tally_comment;
00430 
00431     // Get tally particle
00432     file.getline( line, 100 );
00433     a      = line;
00434     result = get_tally_particle( a, debug, tally_particle );
00435     if( MB_FAILURE == result )
00436     {
00437         // If this line does not specify the particle type, then it is a tally comment.
00438         // Get the comment, then get the particle type from the next line.
00439         tally_comment = line;
00440         file.getline( line, 100 );
00441         a      = line;
00442         result = get_tally_particle( a, debug, tally_particle );
00443         if( MB_SUCCESS != result ) return result;
00444     }
00445     if( debug ) std::cout << "tally_comment=| " << tally_comment << std::endl;
00446 
00447     return MB_SUCCESS;
00448 }
00449 
00450 ErrorCode ReadMCNP5::get_tally_particle( std::string a, bool debug, particle& tally_particle )
00451 {
00452     if( std::string::npos != a.find( "This is a neutron mesh tally." ) ) { tally_particle = NEUTRON; }
00453     else if( std::string::npos != a.find( "This is a photon mesh tally." ) )
00454     {
00455         tally_particle = PHOTON;
00456     }
00457     else if( std::string::npos != a.find( "This is an electron mesh tally." ) )
00458     {
00459         tally_particle = ELECTRON;
00460     }
00461     else
00462         return MB_FAILURE;
00463 
00464     if( debug ) std::cout << "tally_particle=| " << tally_particle << std::endl;
00465 
00466     return MB_SUCCESS;
00467 }
00468 
00469 ErrorCode ReadMCNP5::read_mesh_planes( std::fstream& file, bool debug, std::vector< double > planes[3],
00470                                        coordinate_system& coord_sys )
00471 {
00472     // Tally bin boundaries:
00473     ErrorCode result;
00474     char line[10000];
00475     file.getline( line, 10000 );
00476     std::string a = line;
00477     if( std::string::npos == a.find( "Tally bin boundaries:" ) ) return MB_FAILURE;
00478 
00479     // Decide what coordinate system the tally is using
00480     // first check for Cylindrical coordinates:
00481     file.getline( line, 10000 );
00482     a                        = line;
00483     std::string::size_type b = a.find( "Cylinder origin at" );
00484     if( std::string::npos != b )
00485     {
00486         coord_sys = CYLINDRICAL;
00487         if( debug ) std::cout << "origin, axis, direction=| " << a << std::endl;
00488         std::istringstream ss( a.substr( b + sizeof( "Cylinder origin at" ), 10000 ) );
00489         // The meshtal file does not contain sufficient information to transform
00490         // the particle. Although origin, axs, and vec is needed only origin and
00491         // axs appear in the meshtal file. Why was vec omitted?.
00492         // get origin (not used)
00493         // Cylinder origin at   0.00E+00  0.00E+00  0.00E+00,
00494         // axis in  0.000E+00 0.000E+00 1.000E+00 direction
00495         double origin[3];
00496         if( debug ) std::cout << "origin=| ";
00497         for( int i = 0; i < 3; i++ )
00498         {
00499             ss >> origin[i];
00500             if( debug ) std::cout << origin[i] << " ";
00501         }
00502         if( debug ) std::cout << std::endl;
00503         int length_of_string = 10;
00504         ss.ignore( length_of_string, ' ' );
00505         ss.ignore( length_of_string, ' ' );
00506         ss.ignore( length_of_string, ' ' );
00507         // Get axis (not used)
00508         double axis[3];
00509         if( debug ) std::cout << "axis=| ";
00510         for( int i = 0; i < 3; i++ )
00511         {
00512             ss >> axis[i];
00513             if( debug ) std::cout << axis[i] << " ";
00514         }
00515         if( debug ) std::cout << std::endl;
00516         file.getline( line, 10000 );
00517         a = line;
00518 
00519         // Get r planes
00520         if( debug ) std::cout << "R direction:=| ";
00521         b = a.find( "R direction:" );
00522         if( std::string::npos != b )
00523         {
00524             std::istringstream ss2( a.substr( b + sizeof( "R direction" ), 10000 ) );
00525             result = get_mesh_plane( ss2, debug, planes[0] );
00526             if( MB_SUCCESS != result ) return result;
00527         }
00528         else
00529             return MB_FAILURE;
00530 
00531         // Get z planes
00532         file.getline( line, 10000 );
00533         a = line;
00534         if( debug ) std::cout << "Z direction:=| ";
00535         b = a.find( "Z direction:" );
00536         if( std::string::npos != b )
00537         {
00538             std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
00539             result = get_mesh_plane( ss2, debug, planes[1] );
00540             if( MB_SUCCESS != result ) return result;
00541         }
00542         else
00543             return MB_FAILURE;
00544 
00545         // Get theta planes
00546         file.getline( line, 10000 );
00547         a = line;
00548         if( debug ) std::cout << "Theta direction:=| ";
00549         b = a.find( "Theta direction (revolutions):" );
00550         if( std::string::npos != b )
00551         {
00552             std::istringstream ss2( a.substr( b + sizeof( "Theta direction (revolutions):" ), 10000 ) );
00553             result = get_mesh_plane( ss2, debug, planes[2] );
00554             if( MB_SUCCESS != result ) return result;
00555         }
00556         else
00557             return MB_FAILURE;
00558 
00559         // Cartesian coordinate system:
00560     }
00561     else if( std::string::npos != a.find( "X direction:" ) )
00562     {
00563         coord_sys = CARTESIAN;
00564         // Get x planes
00565         if( debug ) std::cout << "X direction:=| ";
00566         b = a.find( "X direction:" );
00567         if( std::string::npos != b )
00568         {
00569             std::istringstream ss2( a.substr( b + sizeof( "X direction" ), 10000 ) );
00570             result = get_mesh_plane( ss2, debug, planes[0] );
00571             if( MB_SUCCESS != result ) return result;
00572         }
00573         else
00574             return MB_FAILURE;
00575 
00576         // Get y planes
00577         file.getline( line, 10000 );
00578         a = line;
00579         if( debug ) std::cout << "Y direction:=| ";
00580         b = a.find( "Y direction:" );
00581         if( std::string::npos != b )
00582         {
00583             std::istringstream ss2( a.substr( b + sizeof( "Y direction" ), 10000 ) );
00584             result = get_mesh_plane( ss2, debug, planes[1] );
00585             if( MB_SUCCESS != result ) return result;
00586         }
00587         else
00588             return MB_FAILURE;
00589 
00590         // Get z planes
00591         file.getline( line, 10000 );
00592         a = line;
00593         if( debug ) std::cout << "Z direction:=| ";
00594         b = a.find( "Z direction:" );
00595         if( std::string::npos != b )
00596         {
00597             std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
00598             result = get_mesh_plane( ss2, debug, planes[2] );
00599             if( MB_SUCCESS != result ) return result;
00600         }
00601         else
00602             return MB_FAILURE;
00603 
00604         // Spherical coordinate system not yet implemented:
00605     }
00606     else
00607         return MB_FAILURE;
00608 
00609     return MB_SUCCESS;
00610 }
00611 
00612 // Given a stringstream, return a vector of values in the string.
00613 ErrorCode ReadMCNP5::get_mesh_plane( std::istringstream& ss, bool debug, std::vector< double >& plane )
00614 {
00615     double value;
00616     plane.clear();
00617     while( !ss.eof() )
00618     {
00619         ss >> value;
00620         plane.push_back( value );
00621         if( debug ) std::cout << value << " ";
00622     }
00623     if( debug ) std::cout << std::endl;
00624 
00625     return MB_SUCCESS;
00626 }
00627 
00628 ErrorCode ReadMCNP5::read_element_values_and_errors( std::fstream& file, bool /* debug */,
00629                                                      std::vector< double > planes[3], unsigned int n_chopped_x0_planes,
00630                                                      unsigned int n_chopped_x2_planes, particle tally_particle,
00631                                                      double values[], double errors[] )
00632 {
00633     unsigned int index = 0;
00634     // Need to read every line in the file, even if we chop off some elements
00635     for( unsigned int i = 0; i < planes[0].size() - 1 + n_chopped_x0_planes; i++ )
00636     {
00637         for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
00638         {
00639             for( unsigned int k = 0; k < planes[2].size() - 1 + n_chopped_x2_planes; k++ )
00640             {
00641                 char line[100];
00642                 file.getline( line, 100 );
00643                 // If this element has been chopped off, skip it
00644                 if( i < n_chopped_x0_planes ) continue;
00645                 if( k >= planes[2].size() - 1 && k < planes[2].size() - 1 + n_chopped_x2_planes ) continue;
00646                 std::string a = line;
00647                 std::stringstream ss( a );
00648                 double centroid[3];
00649                 double energy;
00650                 // For some reason, photon tallies print the energy in the first column
00651                 if( PHOTON == tally_particle ) ss >> energy;
00652                 // The centroid is not used in this reader
00653                 ss >> centroid[0];
00654                 ss >> centroid[1];
00655                 ss >> centroid[2];
00656                 // Only the tally values and errors are used
00657                 ss >> values[index];
00658                 ss >> errors[index];
00659 
00660                 // Make sure that input data is good
00661                 if( !Util::is_finite( errors[index] ) )
00662                 {
00663                     std::cerr << "found nan error while reading file" << std::endl;
00664                     errors[index] = 1.0;
00665                 }
00666                 if( !Util::is_finite( values[index] ) )
00667                 {
00668                     std::cerr << "found nan value while reading file" << std::endl;
00669                     values[index] = 0.0;
00670                 }
00671 
00672                 index++;
00673             }
00674         }
00675     }
00676 
00677     return MB_SUCCESS;
00678 }
00679 
00680 ErrorCode ReadMCNP5::set_tally_tags( EntityHandle tally_meshset, unsigned int tally_number, char tally_comment[100],
00681                                      particle tally_particle, coordinate_system tally_coord_sys, Tag tally_number_tag,
00682                                      Tag tally_comment_tag, Tag tally_particle_tag, Tag tally_coord_sys_tag )
00683 {
00684     ErrorCode result;
00685     result = MBI->tag_set_data( tally_number_tag, &tally_meshset, 1, &tally_number );
00686     if( MB_SUCCESS != result ) return result;
00687     result = MBI->tag_set_data( tally_comment_tag, &tally_meshset, 1, &tally_comment );
00688     if( MB_SUCCESS != result ) return result;
00689     result = MBI->tag_set_data( tally_particle_tag, &tally_meshset, 1, &tally_particle );
00690     if( MB_SUCCESS != result ) return result;
00691     result = MBI->tag_set_data( tally_coord_sys_tag, &tally_meshset, 1, &tally_coord_sys );
00692     if( MB_SUCCESS != result ) return result;
00693 
00694     return MB_SUCCESS;
00695 }
00696 
00697 ErrorCode ReadMCNP5::create_vertices( std::vector< double > planes[3], bool debug, EntityHandle& start_vert,
00698                                       coordinate_system coord_sys, EntityHandle tally_meshset )
00699 {
00700     // The only info needed to build elements is the mesh plane boundaries.
00701     ErrorCode result;
00702     int n_verts = planes[0].size() * planes[1].size() * planes[2].size();
00703     if( debug ) std::cout << "n_verts=" << n_verts << std::endl;
00704     std::vector< double* > coord_arrays( 3 );
00705     result = readMeshIface->get_node_coords( 3, n_verts, MB_START_ID, start_vert, coord_arrays );
00706     if( MB_SUCCESS != result ) return result;
00707     assert( 0 != start_vert );  // Check for NULL
00708 
00709     for( unsigned int k = 0; k < planes[2].size(); k++ )
00710     {
00711         for( unsigned int j = 0; j < planes[1].size(); j++ )
00712         {
00713             for( unsigned int i = 0; i < planes[0].size(); i++ )
00714             {
00715                 unsigned int idx = ( k * planes[0].size() * planes[1].size() + j * planes[0].size() + i );
00716                 double in[3], out[3];
00717 
00718                 in[0]  = planes[0][i];
00719                 in[1]  = planes[1][j];
00720                 in[2]  = planes[2][k];
00721                 result = transform_point_to_cartesian( in, out, coord_sys );
00722                 if( MB_SUCCESS != result ) return result;
00723 
00724                 // Cppcheck warning (false positive): variable coord_arrays is assigned a value that
00725                 // is never used
00726                 coord_arrays[0][idx] = out[0];
00727                 coord_arrays[1][idx] = out[1];
00728                 coord_arrays[2][idx] = out[2];
00729             }
00730         }
00731     }
00732     Range vert_range( start_vert, start_vert + n_verts - 1 );
00733     result = MBI->add_entities( tally_meshset, vert_range );
00734     if( MB_SUCCESS != result ) return result;
00735 
00736     if( fileIDTag )
00737     {
00738         result = readMeshIface->assign_ids( *fileIDTag, vert_range, nodeId );
00739         if( MB_SUCCESS != result ) return result;
00740         nodeId += vert_range.size();
00741     }
00742 
00743     return MB_SUCCESS;
00744 }
00745 
00746 ErrorCode ReadMCNP5::create_elements( bool debug, std::vector< double > planes[3], unsigned int /*n_chopped_x0_planes*/,
00747                                       unsigned int /*n_chopped_x2_planes*/, EntityHandle start_vert, double* values,
00748                                       double* errors, Tag tally_tag, Tag error_tag, EntityHandle tally_meshset,
00749                                       coordinate_system tally_coord_sys )
00750 {
00751     ErrorCode result;
00752     unsigned int index;
00753     EntityHandle start_element = 0;
00754     unsigned int n_elements    = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
00755     EntityHandle* connect;
00756     result = readMeshIface->get_element_connect( n_elements, 8, MBHEX, MB_START_ID, start_element, connect );
00757     if( MB_SUCCESS != result ) return result;
00758     assert( 0 != start_element );  // Check for NULL
00759 
00760     unsigned int counter = 0;
00761     for( unsigned int i = 0; i < planes[0].size() - 1; i++ )
00762     {
00763         for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
00764         {
00765             for( unsigned int k = 0; k < planes[2].size() - 1; k++ )
00766             {
00767                 index = start_vert + i + j * planes[0].size() + k * planes[0].size() * planes[1].size();
00768                 // For rectangular mesh, the file prints: x y z
00769                 // z changes the fastest and x changes the slowest.
00770                 // This means that the connectivity ordering is not consistent between
00771                 // rectangular and cylindrical mesh.
00772                 if( CARTESIAN == tally_coord_sys )
00773                 {
00774                     connect[0] = index;
00775                     connect[1] = index + 1;
00776                     connect[2] = index + 1 + planes[0].size();
00777                     connect[3] = index + planes[0].size();
00778                     connect[4] = index + planes[0].size() * planes[1].size();
00779                     connect[5] = index + 1 + planes[0].size() * planes[1].size();
00780                     connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
00781                     connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
00782                     // For cylindrical mesh, the file prints: r z theta
00783                     // Theta changes the fastest and r changes the slowest.
00784                 }
00785                 else if( CYLINDRICAL == tally_coord_sys )
00786                 {
00787                     connect[0] = index;
00788                     connect[1] = index + 1;
00789                     connect[2] = index + 1 + planes[0].size() * planes[1].size();
00790                     connect[3] = index + planes[0].size() * planes[1].size();
00791                     connect[4] = index + planes[0].size();
00792                     connect[5] = index + 1 + planes[0].size();
00793                     connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
00794                     connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
00795                 }
00796                 else
00797                     return MB_NOT_IMPLEMENTED;
00798 
00799                 connect += 8;
00800                 counter++;
00801             }
00802         }
00803     }
00804     if( counter != n_elements ) std::cout << "counter=" << counter << " n_elements=" << n_elements << std::endl;
00805 
00806     Range element_range( start_element, start_element + n_elements - 1 );
00807     result = MBI->tag_set_data( tally_tag, element_range, values );
00808     if( MB_SUCCESS != result ) return result;
00809     result = MBI->tag_set_data( error_tag, element_range, errors );
00810     if( MB_SUCCESS != result ) return result;
00811 
00812     // Add the elements to the tally set
00813     result = MBI->add_entities( tally_meshset, element_range );
00814     if( MB_SUCCESS != result ) return result;
00815     if( debug ) std::cout << "Read " << n_elements << " elements from tally." << std::endl;
00816 
00817     if( fileIDTag )
00818     {
00819         result = readMeshIface->assign_ids( *fileIDTag, element_range, elemId );
00820         if( MB_SUCCESS != result ) return result;
00821         elemId += element_range.size();
00822     }
00823 
00824     return MB_SUCCESS;
00825 }
00826 
00827 // Average a tally that was recently read in with one that already exists in
00828 // the interface. Only the existing values will be updated.
00829 ErrorCode ReadMCNP5::average_with_existing_tally( bool debug, unsigned long int& new_nps, unsigned long int nps1,
00830                                                   unsigned int tally_number, Tag tally_number_tag, Tag nps_tag,
00831                                                   Tag tally_tag, Tag error_tag, double* values1, double* errors1,
00832                                                   unsigned int n_elements )
00833 {
00834     // Get the tally number
00835     ErrorCode result;
00836 
00837     // Match the tally number with one from the existing meshtal file
00838     Range matching_tally_number_sets;
00839     const void* const tally_number_val[] = { &tally_number };
00840     result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &tally_number_tag, tally_number_val, 1,
00841                                                 matching_tally_number_sets );
00842     if( MB_SUCCESS != result ) return result;
00843     if( debug ) std::cout << "number of matching meshsets=" << matching_tally_number_sets.size() << std::endl;
00844     assert( 1 == matching_tally_number_sets.size() );
00845 
00846     // Identify which of the meshsets is existing
00847     EntityHandle existing_meshset;
00848     existing_meshset = matching_tally_number_sets.front();
00849 
00850     // Get the existing elements from the set
00851     Range existing_elements;
00852     result = MBI->get_entities_by_type( existing_meshset, MBHEX, existing_elements );
00853     if( MB_SUCCESS != result ) return result;
00854     assert( existing_elements.size() == n_elements );
00855 
00856     // Get the nps of the existing and new tally
00857     unsigned long int nps0;
00858     Range sets_with_this_tag;
00859     result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, sets_with_this_tag );
00860     if( MB_SUCCESS != result ) return result;
00861     if( debug ) std::cout << "number of nps sets=" << sets_with_this_tag.size() << std::endl;
00862     assert( 1 == sets_with_this_tag.size() );
00863     result = MBI->tag_get_data( nps_tag, &sets_with_this_tag.front(), 1, &nps0 );
00864     if( MB_SUCCESS != result ) return result;
00865     if( debug ) std::cout << "nps0=" << nps0 << " nps1=" << nps1 << std::endl;
00866     new_nps = nps0 + nps1;
00867 
00868     // Get tally values from the existing elements
00869     double* values0 = new double[existing_elements.size()];
00870     double* errors0 = new double[existing_elements.size()];
00871     result          = MBI->tag_get_data( tally_tag, existing_elements, values0 );
00872     if( MB_SUCCESS != result )
00873     {
00874         delete[] values0;
00875         delete[] errors0;
00876         return result;
00877     }
00878     result = MBI->tag_get_data( error_tag, existing_elements, errors0 );
00879     if( MB_SUCCESS != result )
00880     {
00881         delete[] values0;
00882         delete[] errors0;
00883         return result;
00884     }
00885 
00886     // Average the values and errors
00887     result = average_tally_values( nps0, nps1, values0, values1, errors0, errors1, n_elements );
00888     if( MB_SUCCESS != result )
00889     {
00890         delete[] values0;
00891         delete[] errors0;
00892         return result;
00893     }
00894 
00895     // Set the averaged information back onto the existing elements
00896     result = MBI->tag_set_data( tally_tag, existing_elements, values0 );
00897     if( MB_SUCCESS != result )
00898     {
00899         delete[] values0;
00900         delete[] errors0;
00901         return result;
00902     }
00903     result = MBI->tag_set_data( error_tag, existing_elements, errors0 );
00904     if( MB_SUCCESS != result )
00905     {
00906         delete[] values0;
00907         delete[] errors0;
00908         return result;
00909     }
00910 
00911     // Cleanup
00912     delete[] values0;
00913     delete[] errors0;
00914 
00915     return MB_SUCCESS;
00916 }
00917 
00918 ErrorCode ReadMCNP5::transform_point_to_cartesian( double* in, double* out, coordinate_system coord_sys )
00919 {
00920     // Transform coordinate system
00921     switch( coord_sys )
00922     {
00923         case CARTESIAN:
00924             out[0] = in[0];  // x
00925             out[1] = in[1];  // y
00926             out[2] = in[2];  // z
00927             break;
00928         // theta is in rotations
00929         case CYLINDRICAL:
00930             out[0] = in[0] * cos( 2 * PI * in[2] );  // x
00931             out[1] = in[0] * sin( 2 * PI * in[2] );  // y
00932             out[2] = in[1];                          // z
00933             break;
00934         case SPHERICAL:
00935             return MB_NOT_IMPLEMENTED;
00936         default:
00937             return MB_NOT_IMPLEMENTED;
00938     }
00939 
00940     return MB_SUCCESS;
00941 }
00942 
00943 // Average two tally values and their error. Return average values in the
00944 // place of first tally values.
00945 ErrorCode ReadMCNP5::average_tally_values( const unsigned long int nps0, const unsigned long int nps1, double* values0,
00946                                            const double* values1, double* errors0, const double* errors1,
00947                                            const unsigned long int n_values )
00948 {
00949     for( unsigned long int i = 0; i < n_values; i++ )
00950     {
00951         // std::cout << " values0=" << values0[i] << " values1=" << values1[i]
00952         //          << " errors0=" << errors0[i] << " errors1=" << errors1[i]
00953         //          << " nps0=" << nps0 << " nps1=" << nps1 << std::endl;
00954         errors0[i] = sqrt( pow( values0[i] * errors0[i] * nps0, 2 ) + pow( values1[i] * errors1[i] * nps1, 2 ) ) /
00955                      ( values0[i] * nps0 + values1[i] * nps1 );
00956 
00957         // It is possible to introduce nans if the values are zero.
00958         if( !Util::is_finite( errors0[i] ) ) errors0[i] = 1.0;
00959 
00960         values0[i] = ( values0[i] * nps0 + values1[i] * nps1 ) / ( nps0 + nps1 );
00961 
00962         // std::cout << " values0=" << values0[i] << " errors0=" << errors0[i] << std::endl;
00963     }
00964     // REMEMBER TO UPDATE NPS0 = NPS0 + NPS1 after this
00965 
00966     return MB_SUCCESS;
00967 }
00968 
00969 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines