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