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