Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
main.cpp File Reference
#include <iostream>
#include <fstream>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <cassert>
#include "mcnpmit.hpp"
#include "moab/CartVect.hpp"
#include "moab/Core.hpp"
#include "MBTagConventions.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/GeomUtil.hpp"
#include "moab/FileOptions.hpp"
#include "../tools/mbcoupler/ElemUtil.hpp"
+ Include dependency graph for main.cpp:

Go to the source code of this file.

Macros

#define MBI   mb_instance()
 

Functions

McnpDatamc_instance ()
 
moab::Interfacemb_instance ()
 
MCNPError read_files (int, char **)
 
MCNPError next_double (std::string, double &, int &)
 
MCNPError next_int (std::string, int &, int &)
 
int main (int argc, char **argv)
 

Variables

moab::Tag coord_tag
 
moab::Tag rotation_tag
 
moab::Tag cfd_heating_tag
 
moab::Tag cfd_error_tag
 
std::string h5m_filename
 
std::string CAD_filename
 
std::string output_filename
 
bool skip_build = false
 
bool read_qnv = false
 
clock_t start_time
 
clock_t load_time
 
clock_t build_time
 
clock_t interp_time
 

Macro Definition Documentation

◆ MBI

#define MBI   mb_instance()

Definition at line 18 of file main.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 37 of file main.cpp.

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,
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,
173  MBresult = MBI->tag_get_handle( "error_tag", 1, moab::MB_TYPE_DOUBLE, cfd_error_tag,
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;
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 }

References moab::Range::begin(), box_max(), box_min(), build_time, moab::AdaptiveKDTree::build_tree(), CAD_filename, CARTESIAN, cfd_error_tag, cfd_heating_tag, moab::Range::clear(), coord_tag, CYLINDRICAL, moab::Range::end(), ErrorCode, moab::Tree::find_all_trees(), h5m_filename, interp_time, load_time, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBHEX, MBI, MBVERTEX, MCNP, MCNP_FAILURE, MCNP_SUCCESS, next_int(), nl, output_filename, moab::ElemUtil::point_in_trilinear_hex(), moab::AdaptiveKDTree::point_search(), read_files(), read_qnv, rotation_tag, moab::Range::size(), skip_build, and start_time.

◆ mb_instance()

moab::Interface * mb_instance ( )

Definition at line 453 of file main.cpp.

454 {
455  static moab::Core inst;
456  return &inst;
457 }

◆ mc_instance()

McnpData * mc_instance ( )

Definition at line 447 of file main.cpp.

448 {
449  static McnpData inst;
450  return &inst;
451 }

◆ next_double()

MCNPError next_double ( std::string  s,
double &  d,
int &  p 
)

Definition at line 399 of file main.cpp.

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 }

References DONE, and MCNP_SUCCESS.

◆ next_int()

MCNPError next_int ( std::string  s,
int &  k,
int &  p 
)

Definition at line 423 of file main.cpp.

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 }

References DONE, and MCNP_SUCCESS.

Referenced by main().

◆ read_files()

MCNPError read_files ( int  argc,
char **  argv 
)

Definition at line 358 of file main.cpp.

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 }

References CAD_filename, h5m_filename, MCNP, MCNP_FAILURE, output_filename, read_qnv, and skip_build.

Referenced by main().

Variable Documentation

◆ build_time

clock_t build_time

Definition at line 35 of file main.cpp.

Referenced by main().

◆ CAD_filename

std::string CAD_filename

Definition at line 29 of file main.cpp.

Referenced by main(), and read_files().

◆ cfd_error_tag

moab::Tag cfd_error_tag

Definition at line 26 of file main.cpp.

Referenced by main().

◆ cfd_heating_tag

moab::Tag cfd_heating_tag

Definition at line 26 of file main.cpp.

Referenced by main().

◆ coord_tag

moab::Tag coord_tag

Definition at line 26 of file main.cpp.

Referenced by main().

◆ h5m_filename

std::string h5m_filename

Definition at line 28 of file main.cpp.

Referenced by main(), and read_files().

◆ interp_time

clock_t interp_time

Definition at line 35 of file main.cpp.

Referenced by main().

◆ load_time

clock_t load_time

Definition at line 35 of file main.cpp.

Referenced by main().

◆ output_filename

std::string output_filename

Definition at line 30 of file main.cpp.

Referenced by main(), and read_files().

◆ read_qnv

bool read_qnv = false

Definition at line 33 of file main.cpp.

Referenced by main(), and read_files().

◆ rotation_tag

moab::Tag rotation_tag

Definition at line 26 of file main.cpp.

Referenced by main().

◆ skip_build

bool skip_build = false

Definition at line 32 of file main.cpp.

Referenced by main(), and read_files().

◆ start_time

clock_t start_time

Definition at line 35 of file main.cpp.

Referenced by main().