 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
40 ReaderIface* ReadMCNP5::factory( Interface* iface )
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 );
49 MBI->query_interface( readMeshIface );
50 assert( NULL != readMeshIface );
51 }
52
53 // Destructor
54 ReadMCNP5::~ReadMCNP5()
55 {
56 if( readMeshIface )
57 {
58 MBI->release_interface( readMeshIface );
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.
131 ErrorCode ReadMCNP5::load_one_file( const char* fname,
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
335 ErrorCode ReadMCNP5::create_tags( Tag& date_and_time_tag,
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,
347 MB_TAG_SPARSE | MB_TAG_CREAT );
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,
352 MB_TAG_SPARSE | MB_TAG_CREAT );
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,
358 MB_TAG_SPARSE | MB_TAG_CREAT );
359 if( MB_SUCCESS != result ) return result;
360 result = MBI->tag_get_handle( "TALLY_PARTICLE_TAG", sizeof( particle ), MB_TYPE_OPAQUE, tally_particle_tag,
361 MB_TAG_SPARSE | MB_TAG_CREAT );
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
374 ErrorCode ReadMCNP5::read_file_header( std::fstream& file,
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
411 ErrorCode ReadMCNP5::set_header_tags( EntityHandle output_meshset,
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
430 ErrorCode ReadMCNP5::read_tally_header( std::fstream& file,
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
501 ErrorCode ReadMCNP5::read_mesh_planes( std::fstream& file,
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
662 ErrorCode ReadMCNP5::read_element_values_and_errors( std::fstream& file,
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
718 ErrorCode ReadMCNP5::set_tally_tags( EntityHandle tally_meshset,
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
793 ErrorCode ReadMCNP5::create_elements( bool debug,
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.
883 ErrorCode ReadMCNP5::average_with_existing_tally( bool debug,
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
979 ErrorCode ReadMCNP5::transform_point_to_cartesian( double* in, double* out, coordinate_system coord_sys )
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
1.9.1.