Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadMCNP5.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #include "ReadMCNP5.hpp"
17 #include "moab/Interface.hpp"
18 #include "moab/ReadUtilIface.hpp"
19 #include "Internals.hpp" // For MB_START_ID
20 #include "moab/Range.hpp"
21 #include "moab/FileOptions.hpp"
22 #include "moab/Util.hpp"
23 
24 #include <iostream>
25 #include <sstream>
26 #include <fstream>
27 #include <vector>
28 #include <cstdlib>
29 #include <cmath>
30 #include <cassert>
31 
32 namespace moab
33 {
34 
35 // Parameters
36 const double ReadMCNP5::PI = 3.141592653589793;
37 const double ReadMCNP5::C2PI = 0.1591549430918954;
38 const double ReadMCNP5::CPI = 0.3183098861837907;
39 
41 {
42  return new ReadMCNP5( iface );
43 }
44 
45 // Constructor
46 ReadMCNP5::ReadMCNP5( Interface* impl ) : MBI( impl ), fileIDTag( NULL ), nodeId( 0 ), elemId( 0 )
47 {
48  assert( NULL != impl );
50  assert( NULL != readMeshIface );
51 }
52 
53 // Destructor
55 {
56  if( readMeshIface )
57  {
59  readMeshIface = 0;
60  }
61 }
62 
63 ErrorCode ReadMCNP5::read_tag_values( const char* /* file_name */,
64  const char* /* tag_name */,
65  const FileOptions& /* opts */,
66  std::vector< int >& /* tag_values_out */,
67  const SubsetList* /* subset_list */ )
68 {
69  return MB_NOT_IMPLEMENTED;
70 }
71 
72 // Load the file as called by the Interface function
73 ErrorCode ReadMCNP5::load_file( const char* filename,
74  const EntityHandle* input_meshset,
75  const FileOptions& options,
76  const ReaderIface::SubsetList* subset_list,
77  const Tag* file_id_tag )
78 {
79  // At this time there is no support for reading a subset of the file
80  if( subset_list )
81  {
82  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for meshtal" );
83  }
84 
85  nodeId = elemId = 0;
86  fileIDTag = file_id_tag;
87 
88  // Average several meshtal files if the AVERAGE_TALLY option is given.
89  // In this case, the integer value is the number of files to average.
90  // If averaging, the filename passed to load_file is assumed to be the
91  // root filename. The files are indexed as "root_filename""index".meshtal.
92  // Indices start with 1.
93  int n_files;
94  bool average = false;
95  ErrorCode result;
96  if( MB_SUCCESS == options.get_int_option( "AVERAGE_TALLY", n_files ) )
97  {
98  // Read the first file (but do not average -> can't average a single file)
99  result = load_one_file( filename, input_meshset, options, average );
100  if( MB_SUCCESS != result ) return result;
101 
102  // Get the root filename
103  std::string root_filename( filename );
104  int length = root_filename.length();
105  root_filename.erase( length - sizeof( ".meshtal" ) );
106 
107  // Average the first file with the rest of the files
108  average = true;
109  for( int i = 2; i <= n_files; i++ )
110  {
111  std::stringstream index;
112  index << i;
113  std::string subsequent_filename = root_filename + index.str() + ".meshtal";
114  result = load_one_file( subsequent_filename.c_str(), input_meshset, options, average );
115  if( MB_SUCCESS != result ) return result;
116  }
117 
118  // If not averaging, read a single file
119  }
120  else
121  {
122  result = load_one_file( filename, input_meshset, options, average );
123  if( MB_SUCCESS != result ) return result;
124  }
125 
126  return MB_SUCCESS;
127 }
128 
129 // This actually reads the file. It creates the mesh elements unless
130 // the file is being averaged with a pre-existing mesh.
132  const EntityHandle* input_meshset,
133  const FileOptions& options,
134  const bool average )
135 {
136  bool debug = false;
137  if( debug ) std::cout << "begin ReadMCNP5::load_one_file" << std::endl;
138 
139  ErrorCode result;
140  std::fstream file;
141  file.open( fname, std::fstream::in );
142  char line[10000];
143 
144  // Create tags
145  Tag date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag, tally_particle_tag,
146  tally_coord_sys_tag, tally_tag, error_tag;
147 
148  result = create_tags( date_and_time_tag, title_tag, nps_tag, tally_number_tag, tally_comment_tag,
149  tally_particle_tag, tally_coord_sys_tag, tally_tag, error_tag );
150  if( MB_SUCCESS != result ) return result;
151 
152  // ******************************************************************
153  // This info exists only at the top of each meshtal file
154  // ******************************************************************
155 
156  // Define characteristics of the entire file
157  char date_and_time[100] = "";
158  char title[100] = "";
159  // This file's number of particles
160  unsigned long int nps;
161  // Sum of this file's and existing file's nps for averaging
162  unsigned long int new_nps;
163 
164  // Read the file header
165  result = read_file_header( file, debug, date_and_time, title, nps );
166  if( MB_SUCCESS != result ) return result;
167 
168  // Blank line
169  file.getline( line, 10000 );
170 
171  // Everything stored in the file being read will be in the input_meshset.
172  // If this is being saved in MOAB, set header tags
173  if( !average && 0 != input_meshset )
174  {
175  result = set_header_tags( *input_meshset, date_and_time, title, nps, date_and_time_tag, title_tag, nps_tag );
176  if( MB_SUCCESS != result ) return result;
177  }
178 
179  // ******************************************************************
180  // This info is repeated for each tally in the meshtal file.
181  // ******************************************************************
182 
183  // If averaging, nps will hold the sum of particles simulated in both tallies.
184  while( !file.eof() )
185  {
186  // Define characteristics of this tally
187  unsigned int tally_number;
188  char tally_comment[100] = "";
189  particle tally_particle;
190  coordinate_system tally_coord_sys;
191  std::vector< double > planes[3];
192  unsigned int n_chopped_x0_planes;
193  unsigned int n_chopped_x2_planes;
194 
195  // Read tally header
196  result = read_tally_header( file, debug, tally_number, tally_comment, tally_particle );
197  if( MB_SUCCESS != result ) return result;
198 
199  // Blank line
200  file.getline( line, 10000 );
201  std::string l = line;
202  // if this string is present then skip the following blank line
203  if( std::string::npos != l.find( "This mesh tally is modified by a dose response function." ) )
204  {
205  file.getline( line, 10000 );
206  }
207 
208  // Read mesh planes
209  result = read_mesh_planes( file, debug, planes, tally_coord_sys );
210  if( MB_SUCCESS != result ) return result;
211 
212  // Get energy boundaries
213  file.getline( line, 10000 );
214  std::string a = line;
215  if( debug ) std::cout << "Energy bin boundaries:=| " << a << std::endl;
216 
217  // Blank
218  file.getline( line, 10000 );
219 
220  // Column headers
221  file.getline( line, 10000 );
222 
223  // If using cylindrical mesh, it may be necessary to chop off the last theta element.
224  // We chop off the last theta plane because the elements will be wrong and skew up
225  // the tree building code. This is
226  // because the hex elements are a linear approximation to the cylindrical elements.
227  // Chopping off the last plane is problem-dependent, and due to MCNP5's mandate
228  // that the cylindrical mesh must span 360 degrees.
229  if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_LAST_AZIMUTHAL_PLANE" ) )
230  {
231  planes[2].pop_back();
232  n_chopped_x2_planes = 1;
233  if( debug ) std::cout << "remove last cylindrical plane option found" << std::endl;
234  }
235  else
236  {
237  n_chopped_x2_planes = 0;
238  }
239 
240  // If using cylindrical mesh, it may be necessary to chop off the first radial element.
241  // These elements extend from the axis and often do not cover areas of interest. For
242  // example, if the mesh was meant to cover r=390-400, the first layer will go from
243  // 0-390 and serve as incorrect source elements for interpolation.
244  if( CYLINDRICAL == tally_coord_sys && MB_SUCCESS == options.get_null_option( "REMOVE_FIRST_RADIAL_PLANE" ) )
245  {
246  std::vector< double >::iterator front = planes[0].begin();
247  planes[0].erase( front );
248  n_chopped_x0_planes = 1;
249  if( debug ) std::cout << "remove first radial plane option found" << std::endl;
250  }
251  else
252  {
253  n_chopped_x0_planes = 0;
254  }
255 
256  // Read the values and errors of each element from the file.
257  // Do not read values that are chopped off.
258  unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
259  double* values = new double[n_elements];
260  double* errors = new double[n_elements];
261  result = read_element_values_and_errors( file, debug, planes, n_chopped_x0_planes, n_chopped_x2_planes,
262  tally_particle, values, errors );
263  if( MB_SUCCESS != result ) return result;
264 
265  // Blank line
266  file.getline( line, 10000 );
267 
268  // ****************************************************************
269  // This tally has been read. If it is not being averaged, build tags,
270  // vertices and elements. If it is being averaged, average the data
271  // with a tally already existing in the MOAB instance.
272  // ****************************************************************
273  if( !average )
274  {
275  EntityHandle tally_meshset;
276  result = MBI->create_meshset( MESHSET_SET, tally_meshset );
277  if( MB_SUCCESS != result ) return result;
278 
279  // Set tags on the tally
280  result = set_tally_tags( tally_meshset, tally_number, tally_comment, tally_particle, tally_coord_sys,
281  tally_number_tag, tally_comment_tag, tally_particle_tag, tally_coord_sys_tag );
282  if( MB_SUCCESS != result ) return result;
283 
284  // The only info needed to build elements is the mesh plane boundaries.
285  // Build vertices...
286  EntityHandle start_vert = 0;
287  result = create_vertices( planes, debug, start_vert, tally_coord_sys, tally_meshset );
288  if( MB_SUCCESS != result ) return result;
289 
290  // Build elements and tag them with tally values and errors, then add
291  // them to the tally_meshset.
292  result = create_elements( debug, planes, n_chopped_x0_planes, n_chopped_x2_planes, start_vert, values,
293  errors, tally_tag, error_tag, tally_meshset, tally_coord_sys );
294  if( MB_SUCCESS != result ) return result;
295 
296  // Add this tally's meshset to the output meshset
297  if( debug ) std::cout << "not averaging tally" << std::endl;
298 
299  // Average the tally values, then delete stuff that was created
300  }
301  else
302  {
303  if( debug ) std::cout << "averaging tally" << std::endl;
304  result = average_with_existing_tally( debug, new_nps, nps, tally_number, tally_number_tag, nps_tag,
305  tally_tag, error_tag, values, errors, n_elements );
306  if( MB_SUCCESS != result ) return result;
307  }
308 
309  // Clean up
310  delete[] values;
311  delete[] errors;
312  }
313 
314  // If we are averaging, delete the remainder of this file's information.
315  // Add the new nps to the existing file's nps if we are averaging.
316  // This is calculated during every tally averaging but only used after the last one.
317  if( average )
318  {
319  Range matching_nps_sets;
320  result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, matching_nps_sets );
321  if( MB_SUCCESS != result ) return result;
322  if( debug ) std::cout << "number of matching nps meshsets=" << matching_nps_sets.size() << std::endl;
323  assert( 1 == matching_nps_sets.size() );
324  result = MBI->tag_set_data( nps_tag, matching_nps_sets, &new_nps );
325  if( MB_SUCCESS != result ) return result;
326 
327  // If this file is not being averaged, return the output meshset.
328  }
329 
330  file.close();
331  return MB_SUCCESS;
332 }
333 
334 // create tags needed for this reader
336  Tag& title_tag,
337  Tag& nps_tag,
338  Tag& tally_number_tag,
339  Tag& tally_comment_tag,
340  Tag& tally_particle_tag,
341  Tag& tally_coord_sys_tag,
342  Tag& tally_tag,
343  Tag& error_tag )
344 {
345  ErrorCode result;
346  result = MBI->tag_get_handle( "DATE_AND_TIME_TAG", 100, MB_TYPE_OPAQUE, date_and_time_tag,
348  if( MB_SUCCESS != result ) return result;
349  result = MBI->tag_get_handle( "TITLE_TAG", 100, MB_TYPE_OPAQUE, title_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
350  if( MB_SUCCESS != result ) return result;
351  result = MBI->tag_get_handle( "NPS_TAG", sizeof( unsigned long int ), MB_TYPE_OPAQUE, nps_tag,
353  if( MB_SUCCESS != result ) return result;
354  result =
355  MBI->tag_get_handle( "TALLY_NUMBER_TAG", 1, MB_TYPE_INTEGER, tally_number_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
356  if( MB_SUCCESS != result ) return result;
357  result = MBI->tag_get_handle( "TALLY_COMMENT_TAG", 100, MB_TYPE_OPAQUE, tally_comment_tag,
359  if( MB_SUCCESS != result ) return result;
360  result = MBI->tag_get_handle( "TALLY_PARTICLE_TAG", sizeof( particle ), MB_TYPE_OPAQUE, tally_particle_tag,
362  if( MB_SUCCESS != result ) return result;
363  result = MBI->tag_get_handle( "TALLY_COORD_SYS_TAG", sizeof( coordinate_system ), MB_TYPE_OPAQUE,
364  tally_coord_sys_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
365  if( MB_SUCCESS != result ) return result;
366  result = MBI->tag_get_handle( "TALLY_TAG", 1, MB_TYPE_DOUBLE, tally_tag, MB_TAG_DENSE | MB_TAG_CREAT );
367  if( MB_SUCCESS != result ) return result;
368  result = MBI->tag_get_handle( "ERROR_TAG", 1, MB_TYPE_DOUBLE, error_tag, MB_TAG_DENSE | MB_TAG_CREAT );
369  if( MB_SUCCESS != result ) return result;
370 
371  return MB_SUCCESS;
372 }
373 
375  bool debug,
376  char date_and_time[100],
377  char title[100],
378  unsigned long int& nps )
379 {
380  // Get simulation date and time
381  // mcnp version 5 ld=11242008 probid = 03/23/09 13:38:56
382  char line[100];
383  file.getline( line, 100 );
384  date_and_time = line;
385  if( debug ) std::cout << "date_and_time=| " << date_and_time << std::endl;
386 
387  // Get simulation title
388  // iter Module 4
389  file.getline( line, 100 );
390  title = line;
391  if( debug ) std::cout << "title=| " << title << std::endl;
392 
393  // Get number of histories
394  // Number of histories used for normalizing tallies = 50000000.00
395  file.getline( line, 100 );
396  std::string a = line;
397  std::string::size_type b = a.find( "Number of histories used for normalizing tallies =" );
398  if( std::string::npos != b )
399  {
400  std::istringstream nps_ss(
401  a.substr( b + sizeof( "Number of histories used for normalizing tallies =" ), 100 ) );
402  nps_ss >> nps;
403  if( debug ) std::cout << "nps=| " << nps << std::endl;
404  }
405  else
406  return MB_FAILURE;
407 
408  return MB_SUCCESS;
409 }
410 
412  char date_and_time[100],
413  char title[100],
414  unsigned long int nps,
415  Tag data_and_time_tag,
416  Tag title_tag,
417  Tag nps_tag )
418 {
419  ErrorCode result;
420  result = MBI->tag_set_data( data_and_time_tag, &output_meshset, 1, &date_and_time );
421  if( MB_SUCCESS != result ) return result;
422  result = MBI->tag_set_data( title_tag, &output_meshset, 1, &title );
423  if( MB_SUCCESS != result ) return result;
424  result = MBI->tag_set_data( nps_tag, &output_meshset, 1, &nps );
425  if( MB_SUCCESS != result ) return result;
426 
427  return MB_SUCCESS;
428 }
429 
431  bool debug,
432  unsigned int& tally_number,
433  char tally_comment[100],
434  particle& tally_particle )
435 {
436  // Get tally number
437  // Mesh Tally Number 104
438  ErrorCode result;
439  char line[100];
440  file.getline( line, 100 );
441  std::string a = line;
442  std::string::size_type b = a.find( "Mesh Tally Number" );
443  if( std::string::npos != b )
444  {
445  std::istringstream tally_number_ss( a.substr( b + sizeof( "Mesh Tally Number" ), 100 ) );
446  tally_number_ss >> tally_number;
447  if( debug ) std::cout << "tally_number=| " << tally_number << std::endl;
448  }
449  else
450  {
451  std::cout << "tally number not found" << std::endl;
452  return MB_FAILURE;
453  }
454 
455  // Next get the tally comment (optional) and particle type
456  // 3mm neutron heating in Be (W/cc)
457  // This is a neutron mesh tally.
458  // std::string tally_comment;
459 
460  // Get tally particle
461  file.getline( line, 100 );
462  a = line;
463  result = get_tally_particle( a, debug, tally_particle );
464  if( MB_FAILURE == result )
465  {
466  // If this line does not specify the particle type, then it is a tally comment.
467  // Get the comment, then get the particle type from the next line.
468  tally_comment = line;
469  file.getline( line, 100 );
470  a = line;
471  result = get_tally_particle( a, debug, tally_particle );
472  if( MB_SUCCESS != result ) return result;
473  }
474  if( debug ) std::cout << "tally_comment=| " << tally_comment << std::endl;
475 
476  return MB_SUCCESS;
477 }
478 
479 ErrorCode ReadMCNP5::get_tally_particle( std::string a, bool debug, particle& tally_particle )
480 {
481  if( std::string::npos != a.find( "This is a neutron mesh tally." ) )
482  {
483  tally_particle = NEUTRON;
484  }
485  else if( std::string::npos != a.find( "This is a photon mesh tally." ) )
486  {
487  tally_particle = PHOTON;
488  }
489  else if( std::string::npos != a.find( "This is an electron mesh tally." ) )
490  {
491  tally_particle = ELECTRON;
492  }
493  else
494  return MB_FAILURE;
495 
496  if( debug ) std::cout << "tally_particle=| " << tally_particle << std::endl;
497 
498  return MB_SUCCESS;
499 }
500 
502  bool debug,
503  std::vector< double > planes[3],
504  coordinate_system& coord_sys )
505 {
506  // Tally bin boundaries:
507  ErrorCode result;
508  char line[10000];
509  file.getline( line, 10000 );
510  std::string a = line;
511  if( std::string::npos == a.find( "Tally bin boundaries:" ) ) return MB_FAILURE;
512 
513  // Decide what coordinate system the tally is using
514  // first check for Cylindrical coordinates:
515  file.getline( line, 10000 );
516  a = line;
517  std::string::size_type b = a.find( "Cylinder origin at" );
518  if( std::string::npos != b )
519  {
520  coord_sys = CYLINDRICAL;
521  if( debug ) std::cout << "origin, axis, direction=| " << a << std::endl;
522  std::istringstream ss( a.substr( b + sizeof( "Cylinder origin at" ), 10000 ) );
523  // The meshtal file does not contain sufficient information to transform
524  // the particle. Although origin, axs, and vec is needed only origin and
525  // axs appear in the meshtal file. Why was vec omitted?.
526  // get origin (not used)
527  // Cylinder origin at 0.00E+00 0.00E+00 0.00E+00,
528  // axis in 0.000E+00 0.000E+00 1.000E+00 direction
529  double origin[3];
530  if( debug ) std::cout << "origin=| ";
531  for( int i = 0; i < 3; i++ )
532  {
533  ss >> origin[i];
534  if( debug ) std::cout << origin[i] << " ";
535  }
536  if( debug ) std::cout << std::endl;
537  int length_of_string = 10;
538  ss.ignore( length_of_string, ' ' );
539  ss.ignore( length_of_string, ' ' );
540  ss.ignore( length_of_string, ' ' );
541  // Get axis (not used)
542  double axis[3];
543  if( debug ) std::cout << "axis=| ";
544  for( int i = 0; i < 3; i++ )
545  {
546  ss >> axis[i];
547  if( debug ) std::cout << axis[i] << " ";
548  }
549  if( debug ) std::cout << std::endl;
550  file.getline( line, 10000 );
551  a = line;
552 
553  // Get r planes
554  if( debug ) std::cout << "R direction:=| ";
555  b = a.find( "R direction:" );
556  if( std::string::npos != b )
557  {
558  std::istringstream ss2( a.substr( b + sizeof( "R direction" ), 10000 ) );
559  result = get_mesh_plane( ss2, debug, planes[0] );
560  if( MB_SUCCESS != result ) return result;
561  }
562  else
563  return MB_FAILURE;
564 
565  // Get z planes
566  file.getline( line, 10000 );
567  a = line;
568  if( debug ) std::cout << "Z direction:=| ";
569  b = a.find( "Z direction:" );
570  if( std::string::npos != b )
571  {
572  std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
573  result = get_mesh_plane( ss2, debug, planes[1] );
574  if( MB_SUCCESS != result ) return result;
575  }
576  else
577  return MB_FAILURE;
578 
579  // Get theta planes
580  file.getline( line, 10000 );
581  a = line;
582  if( debug ) std::cout << "Theta direction:=| ";
583  b = a.find( "Theta direction (revolutions):" );
584  if( std::string::npos != b )
585  {
586  std::istringstream ss2( a.substr( b + sizeof( "Theta direction (revolutions):" ), 10000 ) );
587  result = get_mesh_plane( ss2, debug, planes[2] );
588  if( MB_SUCCESS != result ) return result;
589  }
590  else
591  return MB_FAILURE;
592 
593  // Cartesian coordinate system:
594  }
595  else if( std::string::npos != a.find( "X direction:" ) )
596  {
597  coord_sys = CARTESIAN;
598  // Get x planes
599  if( debug ) std::cout << "X direction:=| ";
600  b = a.find( "X direction:" );
601  if( std::string::npos != b )
602  {
603  std::istringstream ss2( a.substr( b + sizeof( "X direction" ), 10000 ) );
604  result = get_mesh_plane( ss2, debug, planes[0] );
605  if( MB_SUCCESS != result ) return result;
606  }
607  else
608  return MB_FAILURE;
609 
610  // Get y planes
611  file.getline( line, 10000 );
612  a = line;
613  if( debug ) std::cout << "Y direction:=| ";
614  b = a.find( "Y direction:" );
615  if( std::string::npos != b )
616  {
617  std::istringstream ss2( a.substr( b + sizeof( "Y direction" ), 10000 ) );
618  result = get_mesh_plane( ss2, debug, planes[1] );
619  if( MB_SUCCESS != result ) return result;
620  }
621  else
622  return MB_FAILURE;
623 
624  // Get z planes
625  file.getline( line, 10000 );
626  a = line;
627  if( debug ) std::cout << "Z direction:=| ";
628  b = a.find( "Z direction:" );
629  if( std::string::npos != b )
630  {
631  std::istringstream ss2( a.substr( b + sizeof( "Z direction" ), 10000 ) );
632  result = get_mesh_plane( ss2, debug, planes[2] );
633  if( MB_SUCCESS != result ) return result;
634  }
635  else
636  return MB_FAILURE;
637 
638  // Spherical coordinate system not yet implemented:
639  }
640  else
641  return MB_FAILURE;
642 
643  return MB_SUCCESS;
644 }
645 
646 // Given a stringstream, return a vector of values in the string.
647 ErrorCode ReadMCNP5::get_mesh_plane( std::istringstream& ss, bool debug, std::vector< double >& plane )
648 {
649  double value;
650  plane.clear();
651  while( !ss.eof() )
652  {
653  ss >> value;
654  plane.push_back( value );
655  if( debug ) std::cout << value << " ";
656  }
657  if( debug ) std::cout << std::endl;
658 
659  return MB_SUCCESS;
660 }
661 
663  bool /* debug */,
664  std::vector< double > planes[3],
665  unsigned int n_chopped_x0_planes,
666  unsigned int n_chopped_x2_planes,
667  particle tally_particle,
668  double values[],
669  double errors[] )
670 {
671  unsigned int index = 0;
672  // Need to read every line in the file, even if we chop off some elements
673  for( unsigned int i = 0; i < planes[0].size() - 1 + n_chopped_x0_planes; i++ )
674  {
675  for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
676  {
677  for( unsigned int k = 0; k < planes[2].size() - 1 + n_chopped_x2_planes; k++ )
678  {
679  char line[100];
680  file.getline( line, 100 );
681  // If this element has been chopped off, skip it
682  if( i < n_chopped_x0_planes ) continue;
683  if( k >= planes[2].size() - 1 && k < planes[2].size() - 1 + n_chopped_x2_planes ) continue;
684  std::string a = line;
685  std::stringstream ss( a );
686  double centroid[3];
687  double energy;
688  // For some reason, photon tallies print the energy in the first column
689  if( PHOTON == tally_particle ) ss >> energy;
690  // The centroid is not used in this reader
691  ss >> centroid[0];
692  ss >> centroid[1];
693  ss >> centroid[2];
694  // Only the tally values and errors are used
695  ss >> values[index];
696  ss >> errors[index];
697 
698  // Make sure that input data is good
699  if( !Util::is_finite( errors[index] ) )
700  {
701  std::cerr << "found nan error while reading file" << std::endl;
702  errors[index] = 1.0;
703  }
704  if( !Util::is_finite( values[index] ) )
705  {
706  std::cerr << "found nan value while reading file" << std::endl;
707  values[index] = 0.0;
708  }
709 
710  index++;
711  }
712  }
713  }
714 
715  return MB_SUCCESS;
716 }
717 
719  unsigned int tally_number,
720  char tally_comment[100],
721  particle tally_particle,
722  coordinate_system tally_coord_sys,
723  Tag tally_number_tag,
724  Tag tally_comment_tag,
725  Tag tally_particle_tag,
726  Tag tally_coord_sys_tag )
727 {
728  ErrorCode result;
729  result = MBI->tag_set_data( tally_number_tag, &tally_meshset, 1, &tally_number );
730  if( MB_SUCCESS != result ) return result;
731  result = MBI->tag_set_data( tally_comment_tag, &tally_meshset, 1, &tally_comment );
732  if( MB_SUCCESS != result ) return result;
733  result = MBI->tag_set_data( tally_particle_tag, &tally_meshset, 1, &tally_particle );
734  if( MB_SUCCESS != result ) return result;
735  result = MBI->tag_set_data( tally_coord_sys_tag, &tally_meshset, 1, &tally_coord_sys );
736  if( MB_SUCCESS != result ) return result;
737 
738  return MB_SUCCESS;
739 }
740 
741 ErrorCode ReadMCNP5::create_vertices( std::vector< double > planes[3],
742  bool debug,
743  EntityHandle& start_vert,
744  coordinate_system coord_sys,
745  EntityHandle tally_meshset )
746 {
747  // The only info needed to build elements is the mesh plane boundaries.
748  ErrorCode result;
749  int n_verts = planes[0].size() * planes[1].size() * planes[2].size();
750  if( debug ) std::cout << "n_verts=" << n_verts << std::endl;
751  std::vector< double* > coord_arrays( 3 );
752  result = readMeshIface->get_node_coords( 3, n_verts, MB_START_ID, start_vert, coord_arrays );
753  if( MB_SUCCESS != result ) return result;
754  assert( 0 != start_vert ); // Check for NULL
755 
756  for( unsigned int k = 0; k < planes[2].size(); k++ )
757  {
758  for( unsigned int j = 0; j < planes[1].size(); j++ )
759  {
760  for( unsigned int i = 0; i < planes[0].size(); i++ )
761  {
762  unsigned int idx = ( k * planes[0].size() * planes[1].size() + j * planes[0].size() + i );
763  double in[3], out[3];
764 
765  in[0] = planes[0][i];
766  in[1] = planes[1][j];
767  in[2] = planes[2][k];
768  result = transform_point_to_cartesian( in, out, coord_sys );
769  if( MB_SUCCESS != result ) return result;
770 
771  // Cppcheck warning (false positive): variable coord_arrays is assigned a value that
772  // is never used
773  coord_arrays[0][idx] = out[0];
774  coord_arrays[1][idx] = out[1];
775  coord_arrays[2][idx] = out[2];
776  }
777  }
778  }
779  Range vert_range( start_vert, start_vert + n_verts - 1 );
780  result = MBI->add_entities( tally_meshset, vert_range );
781  if( MB_SUCCESS != result ) return result;
782 
783  if( fileIDTag )
784  {
785  result = readMeshIface->assign_ids( *fileIDTag, vert_range, nodeId );
786  if( MB_SUCCESS != result ) return result;
787  nodeId += vert_range.size();
788  }
789 
790  return MB_SUCCESS;
791 }
792 
794  std::vector< double > planes[3],
795  unsigned int /*n_chopped_x0_planes*/,
796  unsigned int /*n_chopped_x2_planes*/,
797  EntityHandle start_vert,
798  double* values,
799  double* errors,
800  Tag tally_tag,
801  Tag error_tag,
802  EntityHandle tally_meshset,
803  coordinate_system tally_coord_sys )
804 {
805  ErrorCode result;
806  unsigned int index;
807  EntityHandle start_element = 0;
808  unsigned int n_elements = ( planes[0].size() - 1 ) * ( planes[1].size() - 1 ) * ( planes[2].size() - 1 );
809  EntityHandle* connect;
810  result = readMeshIface->get_element_connect( n_elements, 8, MBHEX, MB_START_ID, start_element, connect );
811  if( MB_SUCCESS != result ) return result;
812  assert( 0 != start_element ); // Check for NULL
813 
814  unsigned int counter = 0;
815  for( unsigned int i = 0; i < planes[0].size() - 1; i++ )
816  {
817  for( unsigned int j = 0; j < planes[1].size() - 1; j++ )
818  {
819  for( unsigned int k = 0; k < planes[2].size() - 1; k++ )
820  {
821  index = start_vert + i + j * planes[0].size() + k * planes[0].size() * planes[1].size();
822  // For rectangular mesh, the file prints: x y z
823  // z changes the fastest and x changes the slowest.
824  // This means that the connectivity ordering is not consistent between
825  // rectangular and cylindrical mesh.
826  if( CARTESIAN == tally_coord_sys )
827  {
828  connect[0] = index;
829  connect[1] = index + 1;
830  connect[2] = index + 1 + planes[0].size();
831  connect[3] = index + planes[0].size();
832  connect[4] = index + planes[0].size() * planes[1].size();
833  connect[5] = index + 1 + planes[0].size() * planes[1].size();
834  connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
835  connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
836  // For cylindrical mesh, the file prints: r z theta
837  // Theta changes the fastest and r changes the slowest.
838  }
839  else if( CYLINDRICAL == tally_coord_sys )
840  {
841  connect[0] = index;
842  connect[1] = index + 1;
843  connect[2] = index + 1 + planes[0].size() * planes[1].size();
844  connect[3] = index + planes[0].size() * planes[1].size();
845  connect[4] = index + planes[0].size();
846  connect[5] = index + 1 + planes[0].size();
847  connect[6] = index + 1 + planes[0].size() + planes[0].size() * planes[1].size();
848  connect[7] = index + planes[0].size() + planes[0].size() * planes[1].size();
849  }
850  else
851  return MB_NOT_IMPLEMENTED;
852 
853  connect += 8;
854  counter++;
855  }
856  }
857  }
858  if( counter != n_elements ) std::cout << "counter=" << counter << " n_elements=" << n_elements << std::endl;
859 
860  Range element_range( start_element, start_element + n_elements - 1 );
861  result = MBI->tag_set_data( tally_tag, element_range, values );
862  if( MB_SUCCESS != result ) return result;
863  result = MBI->tag_set_data( error_tag, element_range, errors );
864  if( MB_SUCCESS != result ) return result;
865 
866  // Add the elements to the tally set
867  result = MBI->add_entities( tally_meshset, element_range );
868  if( MB_SUCCESS != result ) return result;
869  if( debug ) std::cout << "Read " << n_elements << " elements from tally." << std::endl;
870 
871  if( fileIDTag )
872  {
873  result = readMeshIface->assign_ids( *fileIDTag, element_range, elemId );
874  if( MB_SUCCESS != result ) return result;
875  elemId += element_range.size();
876  }
877 
878  return MB_SUCCESS;
879 }
880 
881 // Average a tally that was recently read in with one that already exists in
882 // the interface. Only the existing values will be updated.
884  unsigned long int& new_nps,
885  unsigned long int nps1,
886  unsigned int tally_number,
887  Tag tally_number_tag,
888  Tag nps_tag,
889  Tag tally_tag,
890  Tag error_tag,
891  double* values1,
892  double* errors1,
893  unsigned int n_elements )
894 {
895  // Get the tally number
896  ErrorCode result;
897 
898  // Match the tally number with one from the existing meshtal file
899  Range matching_tally_number_sets;
900  const void* const tally_number_val[] = { &tally_number };
901  result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &tally_number_tag, tally_number_val, 1,
902  matching_tally_number_sets );
903  if( MB_SUCCESS != result ) return result;
904  if( debug ) std::cout << "number of matching meshsets=" << matching_tally_number_sets.size() << std::endl;
905  assert( 1 == matching_tally_number_sets.size() );
906 
907  // Identify which of the meshsets is existing
908  EntityHandle existing_meshset;
909  existing_meshset = matching_tally_number_sets.front();
910 
911  // Get the existing elements from the set
912  Range existing_elements;
913  result = MBI->get_entities_by_type( existing_meshset, MBHEX, existing_elements );
914  if( MB_SUCCESS != result ) return result;
915  assert( existing_elements.size() == n_elements );
916 
917  // Get the nps of the existing and new tally
918  unsigned long int nps0;
919  Range sets_with_this_tag;
920  result = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &nps_tag, 0, 1, sets_with_this_tag );
921  if( MB_SUCCESS != result ) return result;
922  if( debug ) std::cout << "number of nps sets=" << sets_with_this_tag.size() << std::endl;
923  assert( 1 == sets_with_this_tag.size() );
924  result = MBI->tag_get_data( nps_tag, &sets_with_this_tag.front(), 1, &nps0 );
925  if( MB_SUCCESS != result ) return result;
926  if( debug ) std::cout << "nps0=" << nps0 << " nps1=" << nps1 << std::endl;
927  new_nps = nps0 + nps1;
928 
929  // Get tally values from the existing elements
930  double* values0 = new double[existing_elements.size()];
931  double* errors0 = new double[existing_elements.size()];
932  result = MBI->tag_get_data( tally_tag, existing_elements, values0 );
933  if( MB_SUCCESS != result )
934  {
935  delete[] values0;
936  delete[] errors0;
937  return result;
938  }
939  result = MBI->tag_get_data( error_tag, existing_elements, errors0 );
940  if( MB_SUCCESS != result )
941  {
942  delete[] values0;
943  delete[] errors0;
944  return result;
945  }
946 
947  // Average the values and errors
948  result = average_tally_values( nps0, nps1, values0, values1, errors0, errors1, n_elements );
949  if( MB_SUCCESS != result )
950  {
951  delete[] values0;
952  delete[] errors0;
953  return result;
954  }
955 
956  // Set the averaged information back onto the existing elements
957  result = MBI->tag_set_data( tally_tag, existing_elements, values0 );
958  if( MB_SUCCESS != result )
959  {
960  delete[] values0;
961  delete[] errors0;
962  return result;
963  }
964  result = MBI->tag_set_data( error_tag, existing_elements, errors0 );
965  if( MB_SUCCESS != result )
966  {
967  delete[] values0;
968  delete[] errors0;
969  return result;
970  }
971 
972  // Cleanup
973  delete[] values0;
974  delete[] errors0;
975 
976  return MB_SUCCESS;
977 }
978 
980 {
981  // Transform coordinate system
982  switch( coord_sys )
983  {
984  case CARTESIAN:
985  out[0] = in[0]; // x
986  out[1] = in[1]; // y
987  out[2] = in[2]; // z
988  break;
989  // theta is in rotations
990  case CYLINDRICAL:
991  out[0] = in[0] * cos( 2 * PI * in[2] ); // x
992  out[1] = in[0] * sin( 2 * PI * in[2] ); // y
993  out[2] = in[1]; // z
994  break;
995  case SPHERICAL:
996  return MB_NOT_IMPLEMENTED;
997  default:
998  return MB_NOT_IMPLEMENTED;
999  }
1000 
1001  return MB_SUCCESS;
1002 }
1003 
1004 // Average two tally values and their error. Return average values in the
1005 // place of first tally values.
1006 ErrorCode ReadMCNP5::average_tally_values( const unsigned long int nps0,
1007  const unsigned long int nps1,
1008  double* values0,
1009  const double* values1,
1010  double* errors0,
1011  const double* errors1,
1012  const unsigned long int n_values )
1013 {
1014  for( unsigned long int i = 0; i < n_values; i++ )
1015  {
1016  // std::cout << " values0=" << values0[i] << " values1=" << values1[i]
1017  // << " errors0=" << errors0[i] << " errors1=" << errors1[i]
1018  // << " nps0=" << nps0 << " nps1=" << nps1 << std::endl;
1019  errors0[i] = sqrt( pow( values0[i] * errors0[i] * nps0, 2 ) + pow( values1[i] * errors1[i] * nps1, 2 ) ) /
1020  ( values0[i] * nps0 + values1[i] * nps1 );
1021 
1022  // It is possible to introduce nans if the values are zero.
1023  if( !Util::is_finite( errors0[i] ) ) errors0[i] = 1.0;
1024 
1025  values0[i] = ( values0[i] * nps0 + values1[i] * nps1 ) / ( nps0 + nps1 );
1026 
1027  // std::cout << " values0=" << values0[i] << " errors0=" << errors0[i] << std::endl;
1028  }
1029  // REMEMBER TO UPDATE NPS0 = NPS0 + NPS1 after this
1030 
1031  return MB_SUCCESS;
1032 }
1033 
1034 } // namespace moab