 |
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 #include <iostream>
2 #include <fstream>
3 #include <cmath>
4 #include <ctime>
5 #include <cstdlib>
6 #include <cassert>
7
8 #include "mcnpmit.hpp"
9
10 #include "moab/CartVect.hpp"
11 #include "moab/Core.hpp"
12 #include "MBTagConventions.hpp"
13 #include "moab/AdaptiveKDTree.hpp"
14 #include "moab/GeomUtil.hpp"
15 #include "moab/FileOptions.hpp"
16 #include "../tools/mbcoupler/ElemUtil.hpp"
17
18 #define MBI mb_instance()
19
20 McnpData* mc_instance();
21 moab::Interface* mb_instance();
22 MCNPError read_files( int, char** );
23 MCNPError next_double( std::string, double&, int& );
24 MCNPError next_int( std::string, int&, int& );
25
26 moab::Tag coord_tag, rotation_tag, cfd_heating_tag, cfd_error_tag;
27
28 std::string h5m_filename;
29 std::string CAD_filename;
30 std::string output_filename;
31
32 bool skip_build = false;
33 bool read_qnv = false;
34
35 clock_t start_time, load_time, build_time, interp_time;
36
37 int main( int argc, char** argv )
38 {
39 MCNPError result;
40
41 start_time = clock();
42
43 // Read in file names from command line
44 result = read_files( argc, argv );
45 if( result == MCNP_FAILURE ) return 1;
46
47 result = MCNP->initialize_tags();
48
49 // Parse the MCNP input file
50 if( !skip_build )
51 {
52
53 result = MCNP->read_mcnpfile( skip_build );
54 if( result == MCNP_FAILURE )
55 {
56 std::cout << "Failure reading MCNP file!" << std::endl;
57 return 1;
58 }
59 }
60
61 load_time = clock() - start_time;
62
63 // Make the KD-Tree
64 moab::ErrorCode MBresult;
65 moab::AdaptiveKDTree kdtree( MBI );
66 moab::EntityHandle root;
67
68 MBI->tag_get_handle( "CoordTag", 1, moab::MB_TYPE_INTEGER, coord_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
69 MBI->tag_get_handle( "RotationTag", 16, moab::MB_TYPE_DOUBLE, rotation_tag,
70 moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
71
72 if( skip_build )
73 {
74 MBresult = MBI->load_mesh( h5m_filename.c_str() );
75
76 if( moab::MB_SUCCESS == MBresult )
77 {
78 std::cout << std::endl << "Read in mesh from h5m file." << std::endl << std::endl;
79 std::cout << "Querying mesh file..." << std::endl;
80 }
81 else
82 {
83 std::cout << "Failure reading h5m file!" << std::endl;
84 std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
85 std::string message;
86 if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
87 std::cerr << "Error message: " << message << std::endl;
88 return 1;
89 }
90
91 moab::Range tmprange;
92 kdtree.find_all_trees( tmprange );
93 root = tmprange[0];
94 }
95 else
96 {
97 std::cout << "Building KD-Tree..." << std::endl;
98 moab::FileOptions opts( "CANDIDATE_PLANE_SET=SUBDIVISION" );
99 MBresult = kdtree.build_tree( MCNP->elem_handles, &root, &opts );
100 if( MBresult == moab::MB_SUCCESS )
101 {
102
103 MBI->tag_set_data( coord_tag, &root, 1, &( MCNP->coord_system ) );
104 MBI->tag_set_data( rotation_tag, &root, 1, &( MCNP->rotation_matrix ) );
105
106 std::cout << "KD-Tree has been built successfully!" << std::endl << std::endl;
107 MBresult = MBI->write_mesh( ( MCNP->MCNP_filename + ".h5m" ).c_str() );
108
109 std::cout << "Querying mesh file..." << std::endl;
110 }
111 else
112 {
113 std::cout << "Error building KD-Tree!" << std::endl << std::endl;
114 std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
115 std::string message;
116 if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
117 std::cerr << "Error message: " << message << std::endl;
118 return 1;
119 }
120 }
121
122 int coord_sys;
123 double rmatrix[16];
124
125 MBresult = MBI->tag_get_data( coord_tag, &root, 1, &coord_sys );MB_CHK_ERR( MBresult );
126 MBresult = MBI->tag_get_data( rotation_tag, &root, 1, &rmatrix );MB_CHK_ERR( MBresult );
127
128 build_time = clock() - load_time;
129
130 // Read the CAD mesh data and query the tree
131 std::ifstream cadfile;
132 std::ofstream outfile;
133
134 outfile.open( output_filename.c_str() );
135
136 int num_pts;
137 int n;
138 long int nl = 0;
139 char* ctmp;
140 int elems_read = 0;
141 int p = 0;
142 char line[10000];
143
144 // Used only when reading a mesh file to get vertex info
145 double* cfd_coords = NULL;
146 moab::Range::iterator cfd_iter;
147 moab::EntityHandle meshset;
148
149 if( read_qnv )
150 {
151 cadfile.open( CAD_filename.c_str() );
152 cadfile.getline( line, 10000 );
153 cadfile.getline( line, 10000 );
154 result = next_int( line, num_pts, p );
155 }
156 else
157 {
158
159 meshset = 0;
160 MBresult = MBI->load_file( CAD_filename.c_str(), &meshset );MB_CHK_ERR( MBresult );
161 assert( 0 != meshset );
162
163 moab::Range cfd_verts;
164 MBresult = MBI->get_entities_by_type( meshset, moab::MBVERTEX, cfd_verts, true );MB_CHK_ERR( MBresult );
165 num_pts = cfd_verts.size();
166
167 cfd_coords = new double[3 * num_pts];
168 MBresult = MBI->get_coords( cfd_verts, cfd_coords );MB_CHK_ERR( MBresult );
169
170 cfd_iter = cfd_verts.begin();
171 MBresult = MBI->tag_get_handle( "heating_tag", 1, moab::MB_TYPE_DOUBLE, cfd_heating_tag,
172 moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );
173 MBresult = MBI->tag_get_handle( "error_tag", 1, moab::MB_TYPE_DOUBLE, cfd_error_tag,
174 moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );
175
176 std::cout << std::endl << "Read in mesh with query points." << std::endl << std::endl;
177 }
178
179 double testpt[3];
180 double transformed_pt[3];
181 double taldata;
182 double errdata;
183
184 moab::CartVect testvc;
185
186 bool found = false;
187
188 // MBRange verts;
189 std::vector< moab::EntityHandle > verts;
190 moab::Range range;
191 moab::CartVect box_max, box_min;
192
193 moab::CartVect hexverts[8];
194 moab::CartVect tmp_cartvect;
195 std::vector< double > coords;
196
197 double tal_sum = 0.0, err_sum = 0.0, tal_sum_sqr = 0.0, err_sum_sqr = 0.0;
198
199 // double davg = 0.0;
200 // unsigned int nmax = 0, nmin = 1000000000 ;
201
202 for( unsigned int i = 0; i < (unsigned int)num_pts; i++ )
203 {
204
205 // if (i%status_freq == 0)
206 // std::cerr << "Completed " << i/status_freq << "%" << std::endl;
207
208 // Grab the coordinates to query
209 if( read_qnv )
210 {
211 cadfile.getline( line, 10000 );
212
213 nl = std::strtol( line, &ctmp, 10 );
214 n = (unsigned int)nl;
215 testpt[0] = std::strtod( ctmp + 1, &ctmp );
216 testpt[1] = std::strtod( ctmp + 1, &ctmp );
217 testpt[2] = std::strtod( ctmp + 1, NULL );
218 }
219 else
220 {
221 testpt[0] = cfd_coords[3 * i];
222 testpt[1] = cfd_coords[3 * i + 1];
223 testpt[2] = cfd_coords[3 * i + 2];
224 n = i + 1;
225 }
226
227 result = MCNP->transform_point( testpt, transformed_pt, coord_sys, rmatrix );
228
229 testvc[0] = transformed_pt[0];
230 testvc[1] = transformed_pt[1];
231 testvc[2] = transformed_pt[2];
232
233 // Find the leaf containing the point
234 moab::EntityHandle tree_node;
235 MBresult = kdtree.point_search( transformed_pt, tree_node );
236 if( moab::MB_SUCCESS != MBresult )
237 {
238 double x = 0.0, y = 0.0, z = 0.0;
239 if( CARTESIAN == coord_sys )
240 {
241 x = testvc[0];
242 y = testvc[1];
243 z = testvc[2];
244 }
245 else if( CYLINDRICAL == coord_sys )
246 {
247 x = testvc[0] * cos( 2 * M_PI * testvc[2] );
248 y = testvc[0] * sin( 2 * M_PI * testvc[2] );
249 z = testvc[1];
250 }
251 else
252 {
253 std::cout << "MOAB WARNING: Unhandled error code during point search in KdTree, "
254 "ErrorCode = "
255 << MBresult << " and Coord xyz=" << x << " " << y << " " << z << std::endl;
256 }
257 std::cout << "No leaf found, MCNP coord xyz=" << x << " " << y << " " << z << std::endl;
258 ++cfd_iter;
259 continue;
260 }
261
262 range.clear();
263 MBresult = MBI->get_entities_by_type( tree_node, moab::MBHEX, range );MB_CHK_ERR( MBresult );
264
265 // davg += (double) range.size();
266 // if (range.size() > nmax) nmax = range.size();
267 // if (range.size() < nmin) nmin = range.size();
268
269 for( moab::Range::iterator rit = range.begin(); rit != range.end(); ++rit )
270 {
271 verts.clear();
272 const moab::EntityHandle* connect;
273 int num_connect;
274 MBresult = MBI->get_connectivity( *rit, connect, num_connect, true );MB_CHK_ERR( MBresult );
275
276 coords.resize( 3 * num_connect );
277 MBresult = MBI->get_coords( connect, num_connect, &coords[0] );MB_CHK_ERR( MBresult );
278
279 for( unsigned int j = 0; j < (unsigned int)num_connect; j++ )
280 {
281 hexverts[j][0] = coords[3 * j];
282 hexverts[j][1] = coords[3 * j + 1];
283 hexverts[j][2] = coords[3 * j + 2];
284 }
285
286 if( moab::ElemUtil::point_in_trilinear_hex( hexverts, testvc, 1.e-6 ) )
287 {
288 MBresult = MBI->tag_get_data( MCNP->tally_tag, &( *rit ), 1, &taldata );MB_CHK_ERR( MBresult );
289 MBresult = MBI->tag_get_data( MCNP->relerr_tag, &( *rit ), 1, &errdata );MB_CHK_ERR( MBresult );
290
291 outfile << n << "," << testpt[0] << "," << testpt[1] << "," << testpt[2] << "," << taldata << ","
292 << errdata << std::endl;
293
294 if( !read_qnv )
295 {
296 MBresult = MBI->tag_set_data( cfd_heating_tag, &( *cfd_iter ), 1, &taldata );MB_CHK_ERR( MBresult );
297 MBresult = MBI->tag_set_data( cfd_error_tag, &( *cfd_iter ), 1, &errdata );MB_CHK_ERR( MBresult );
298 }
299
300 found = true;
301 elems_read++;
302
303 tal_sum = tal_sum + taldata;
304 err_sum = err_sum + errdata;
305 tal_sum_sqr = tal_sum_sqr + taldata * taldata;
306 err_sum_sqr = err_sum_sqr + errdata * errdata;
307
308 break;
309 }
310 }
311
312 if( !read_qnv ) ++cfd_iter;
313
314 if( !found )
315 {
316 std::cout << n << " " << testvc << std::endl;
317 }
318 found = false;
319 }
320
321 cadfile.close();
322 outfile.close();
323
324 if( result == MCNP_SUCCESS )
325 {
326 std::cout << "Success! " << elems_read << " elements interpolated." << std::endl << std::endl;
327 }
328 else
329 {
330 std::cout << "Failure during query! " << elems_read << " elements interpolated." << std::endl << std::endl;
331 }
332
333 double tal_std_dev = sqrt( ( 1.0 / elems_read ) * ( tal_sum_sqr - ( 1.0 / elems_read ) * tal_sum * tal_sum ) );
334 double err_std_dev = sqrt( ( 1.0 / elems_read ) * ( err_sum_sqr - ( 1.0 / elems_read ) * err_sum * err_sum ) );
335
336 std::cout << "Tally Mean: " << tal_sum / elems_read << std::endl;
337 std::cout << "Tally Standard Deviation: " << tal_std_dev << std::endl;
338 std::cout << "Error Mean: " << err_sum / elems_read << std::endl;
339 std::cout << "Error Standard Deviation: " << err_std_dev << std::endl;
340
341 interp_time = clock() - build_time;
342
343 if( !read_qnv )
344 {
345 std::string out_mesh_fname = output_filename;
346 MBresult = MBI->write_mesh( ( out_mesh_fname + ".h5m" ).c_str(), &meshset, 1 );
347 // MBresult = MBI->write_file( (cfd_mesh_fname + ".vtk").c_str(), "vtk", NULL, &meshset, 1,
348 // &cfd_heating_tag, 1);
349 }
350
351 std::cout << "Time to read in file: " << (double)load_time / CLOCKS_PER_SEC << std::endl;
352 std::cout << "Time to build kd-tree: " << (double)build_time / CLOCKS_PER_SEC << std::endl;
353 std::cout << "Time to interpolate data: " << (double)interp_time / CLOCKS_PER_SEC << std::endl;
354
355 return 0;
356 }
357
358 MCNPError read_files( int argc, char** argv )
359 {
360 MCNPError result = MCNP_FAILURE;
361
362 // Check to see if appropriate command lines specified
363 if( argc < 3 )
364 {
365 std::cout << "Source and Target mesh filenames NOT specified!";
366 std::cout << std::endl;
367 return MCNP_FAILURE;
368 }
369
370 // Set the MCNP or H5M filename
371 std::string str;
372 str = argv[1];
373
374 unsigned int itmp = str.find( ".h5m" );
375 if( ( itmp > 0 ) && ( itmp < str.length() ) )
376 {
377 skip_build = true;
378 h5m_filename = str;
379 }
380 else
381 {
382 result = MCNP->set_filename( str );
383 }
384
385 // Set the CAD filename
386 str = argv[2];
387 CAD_filename = str;
388
389 itmp = str.find( ".qnv" );
390 if( ( itmp > 0 ) && ( itmp < str.length() ) ) read_qnv = true;
391
392 // Set the output filename
393 str = argv[3];
394 output_filename = str;
395
396 return result;
397 }
398
399 MCNPError next_double( std::string s, double& d, int& p )
400 {
401
402 unsigned int slen = s.length();
403 unsigned int j;
404
405 for( unsigned int i = p; i < slen; i++ )
406 {
407 if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
408 {
409
410 j = s.find( ",", i );
411 if( j > slen ) j = slen;
412
413 d = std::atof( s.substr( i, j - i ).c_str() );
414 p = j + 1;
415
416 return MCNP_SUCCESS;
417 }
418 }
419
420 return DONE;
421 }
422
423 MCNPError next_int( std::string s, int& k, int& p )
424 {
425
426 unsigned int slen = s.length();
427 unsigned int j;
428
429 for( unsigned int i = p; i < slen; i++ )
430 {
431 if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
432 {
433
434 j = s.find( ",", i );
435 if( j > slen ) j = slen;
436
437 k = std::atoi( s.substr( i, j - i ).c_str() );
438 p = j + 1;
439
440 return MCNP_SUCCESS;
441 }
442 }
443
444 return DONE;
445 }
446
447 McnpData* mc_instance()
448 {
449 static McnpData inst;
450 return &inst;
451 }
452
453 moab::Interface* mb_instance()
454 {
455 static moab::Core inst;
456 return &inst;
457 }
1.9.1.