MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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