Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
mcnpmit.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <cstdlib>
4 #include "mcnpmit.hpp"
5 #include "moab/CartVect.hpp"
6 #include <cmath>
7 
9 
10 // Parameters
11 // const double pi = 3.141592653589793;
12 const double c2pi = 0.1591549430918954;
13 // const double cpi = 0.3183098861837907;
14 
15 MCNPError next_number( std::string, double&, int& );
16 int how_many_numbers( std::string );
17 MCNPError read_numbers( std::string, int, std::vector< double >& );
18 
19 // Constructor
21 {
22 
23  // Default value for coordinate system
24  coord_system = 0;
25 
26  // Default rotation matrix is identity matrix
27  for( int i = 0; i < 4; i++ )
28  {
29  for( int j = 0; j < 4; j++ )
30  {
31  if( i == j )
32  rotation_matrix[4 * i + j] = 1;
33  else
34  rotation_matrix[4 * i + j] = 0;
35  }
36  }
37 }
38 
39 // Destructor
41 {
42 
43  // Vertices and elements
44  MCNP_vertices.clear();
45 }
46 
47 // Setting and retrieving coordinate sysem
49 {
50  coord_system = k;
51  return MCNP_SUCCESS;
52 }
54 {
55  return coord_system;
56 }
57 
58 // Setting and retrieving roation matrix
60 {
61  for( int i = 0; i < 16; i++ )
62  {
63  rotation_matrix[i] = r[i];
64  }
65  return MCNP_SUCCESS;
66 }
68 {
69  return rotation_matrix;
70 }
71 
72 // Set the filename
73 MCNPError McnpData::set_filename( std::string fname )
74 {
75  MCNP_filename = fname;
76  return MCNP_SUCCESS;
77 }
79 {
80  return MCNP_filename;
81 }
82 
83 // Reading the MCNP file
85 {
86 
87  MCNPError result;
88  moab::ErrorCode MBresult;
89  moab::CartVect tvect;
90 
91  std::vector< double > xvec[3];
92 
93  // Open the file
94  std::ifstream mcnpfile;
95  mcnpfile.open( MCNP_filename.c_str() );
96  if( !mcnpfile )
97  {
98  std::cout << "Unable to open MCNP data file." << std::endl;
99  return MCNP_FAILURE;
100  }
101  std::cout << std::endl;
102  std::cout << "Reading MCNP input file..." << std::endl;
103 
104  // Prepare for file reading ...
105  char line[10000];
106  int mode = 0; // Set the file reading mode to read proper data
107  int nv[3];
108 
109  // Read in the file ...
110  while( !mcnpfile.eof() )
111  {
112 
113  mcnpfile.getline( line, 10000 );
114  // std::cout << line << std::endl;
115 
116  switch( mode )
117  {
118  case 0: // First line is a title
119  mode++;
120  break;
121  case 1: // Coordinate system
122  mode++;
123  result = read_coord_system( line );
124  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
125  break;
126  case 2: // Rotation matrix
127  mode++;
128  for( int i = 0; i < 4; i++ )
129  {
130  mcnpfile.getline( line, 10000 );
131  result = read_rotation_matrix( line, i );
132  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
133  }
134  if( skip_mesh ) return MCNP_SUCCESS;
135  break;
136  case 3: // Read in vertices and build elements
137  mode++;
138 
139  for( int i = 0; i < 3; i++ )
140  {
141  // How many points in the x[i]-direction
142  nv[i] = how_many_numbers( line );
143  if( nv[i] <= 0 ) return MCNP_FAILURE;
144 
145  // Get space and read in these points
146  result = read_numbers( line, nv[i], xvec[i] );
147  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
148 
149  // Update to the next line
150  mcnpfile.getline( line, 10000 );
151  }
152 
153  // Make the elements and vertices
154  result = make_elements( xvec, nv );
155  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
156  break;
157  case 4: // Read in tally data, make, and tag elements
158  mode++;
159  moab::EntityHandle elemhandle;
160 
161  moab::EntityHandle vstart, vijk;
162  moab::EntityHandle connect[8];
163  // double d[3];
164 
165  // vstart = MCNP_vertices.front();
166  vstart = *( vert_handles.begin() );
167 
168  for( int i = 0; i < nv[0] - 1; i++ )
169  {
170  for( int j = 0; j < nv[1] - 1; j++ )
171  {
172  for( int k = 0; k < nv[2] - 1; k++ )
173  {
174  vijk = vstart + ( i + j * nv[0] + k * nv[0] * nv[1] );
175 
176  // std::cout << vijk << std::endl;
177 
178  connect[0] = vijk;
179  connect[1] = vijk + 1;
180  connect[2] = vijk + 1 + nv[0];
181  connect[3] = vijk + nv[0];
182  connect[4] = vijk + nv[0] * nv[1];
183  connect[5] = vijk + 1 + nv[0] * nv[1];
184  connect[6] = vijk + 1 + nv[0] + nv[0] * nv[1];
185  connect[7] = vijk + nv[0] + nv[0] * nv[1];
186 
187  MBresult = MBI->create_element( moab::MBHEX, connect, 8, elemhandle );
188  if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
189  elem_handles.insert( elemhandle );
190 
191  mcnpfile.getline( line, 10000 );
192  result = extract_tally_data( line, elemhandle );
193  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
194  }
195  }
196  }
197 
198  /*
199  for (MBRange::iterator rit=vert_handles.begin(); rit !=
200  vert_handles.end(); ++rit) { std::cout << *rit << std::endl;
201  }
202 
203 
204  for (int i=0; i < nv[0]-1; i++) {
205  for (int j=0; j < nv[1]-1; j++) {
206  for (int k=0; k < nv[2]-1; k++) {
207  vijk = vstart + (i + j*nv[0] + k*nv[0]*nv[1]);
208  connect[0] = vijk;
209  connect[1] = vijk + 1;
210  connect[2] = vijk + 1 + nv[0];
211  connect[3] = vijk + nv[0];
212  connect[4] = vijk + nv[0]*nv[1];
213  connect[5] = vijk + 1 + nv[0]*nv[1];
214  connect[6] = vijk + 1 + nv[0] + nv[0]*nv[1];
215  connect[7] = vijk + nv[0] + nv[0]*nv[1];
216 
217  MBresult = MBI->create_element(MBHEX, connect, 8,
218  elemhandle); if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
219  elem_handles.insert(elemhandle);
220 
221  mcnpfile.getline(line, 10000);
222  result = extract_tally_data(line, elemhandle);
223  if (result == MCNP_FAILURE) return MCNP_FAILURE;
224 
225  }
226  }
227  }
228  */
229  break;
230  case 5: // Ckeck for weirdness at end of file
231  if( !mcnpfile.eof() ) return MCNP_FAILURE;
232  break;
233  }
234  }
235 
236  std::cout << "SUCCESS! Read in " << elem_handles.size() << " elements!" << std::endl << std::endl;
237  // MCNP_vertices.clear();
239  MCNP_elems.clear();
240  return MCNP_SUCCESS;
241 }
242 
244 {
245 
246  if( ( s.find( "Box" ) < 100 ) || ( s.find( "xyz" ) < 100 ) )
248  else if( s.find( "Cyl" ) < 100 )
250  else if( s.find( "Sph" ) < 100 )
252  else
253  return MCNP_FAILURE;
254 
255  return MCNP_SUCCESS;
256 }
257 
259 {
260 
261  int fpos = 0;
262  MCNPError result;
263 
264  for( int j = 0; j < 4; j++ )
265  {
266  result = next_number( s, rotation_matrix[4 * i + j], fpos );
267  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
268  }
269 
270  return MCNP_SUCCESS;
271 }
272 
273 MCNPError McnpData::make_elements( std::vector< double > x[3], int* n )
274 {
275 
276  // double v[3];
277  // MBEntityHandle dumhandle;
278  // MBEntityHandle vstart, vijk;
279  unsigned int num_verts = n[0] * n[1] * n[2];
280  double* coords;
281  coords = new double[3 * num_verts];
282 
283  /*
284  // Enter the vertices ...
285  for (int k=0; k < n[2]; k++) {
286  v[2] = x[2].at(k);
287  for (int j=0; j < n[1]; j++) {
288  v[1] = x[1].at(j);
289  for (int i=0; i < n[0]; i++) {
290  v[0] = x[0].at(i);
291  MBresult = MBI->create_vertex(v, dumhandle);
292  if (MBresult != MB_SUCCESS) return MCNP_FAILURE;
293  MCNP_vertices.push_back(dumhandle);
294 
295  }
296  }
297  }
298  */
299 
300  // Enter the vertices ...
301  for( int k = 0; k < n[2]; k++ )
302  {
303  for( int j = 0; j < n[1]; j++ )
304  {
305  for( int i = 0; i < n[0]; i++ )
306  {
307  unsigned int ijk = 3 * ( k * n[0] * n[1] + j * n[0] + i );
308  coords[ijk] = x[0][i];
309  coords[ijk + 1] = x[1][j];
310  coords[ijk + 2] = x[2][k];
311 
312  // std::cout << coords[ijk] << " " << coords[ijk+1] << " "
313  // << coords[ijk+2] << std::endl;
314 
315  // MCNP_vertices.push_back(dumhandle);
316  }
317  }
318  }
319 
320  MBI->create_vertices( coords, num_verts, vert_handles );
321 
322  delete[] coords;
323  return MCNP_SUCCESS;
324 }
325 
327 {
328 
331 
332  return MCNP_SUCCESS;
333 }
334 
336 {
337 
338  int fpos = 0;
339  double d = 0;
340 
341  MCNPError result;
342  moab::ErrorCode MBresult;
343 
344  // Discard first three lines
345  for( int i = 0; i < 3; i++ )
346  {
347  result = next_number( s, d, fpos );
348  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
349  }
350  // Need to read in tally entry and tag ...
351  result = next_number( s, d, fpos );
352  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
353  MBresult = MBI->tag_set_data( tally_tag, &handle, 1, &d );
354  if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
355 
356  // Need to read in relative error entry and tag ...
357  result = next_number( s, d, fpos );
358  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
359  MBresult = MBI->tag_set_data( relerr_tag, &handle, 1, &d );
360  if( MBresult != moab::MB_SUCCESS ) return MCNP_FAILURE;
361 
362  return MCNP_SUCCESS;
363 }
364 
365 MCNPError next_number( std::string s, double& d, int& p )
366 {
367 
368  unsigned int slen = s.length();
369  unsigned int j;
370 
371  for( unsigned int i = p; i < slen; i++ )
372  {
373  if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
374  {
375 
376  j = s.find( " ", i );
377 
378  if( j > slen ) j = slen;
379 
380  // Extract the number out of the string
381  d = std::atof( s.substr( i, j - i ).c_str() );
382  p = j + 1;
383 
384  return MCNP_SUCCESS;
385  }
386  }
387 
388  return DONE;
389 }
390 
391 int how_many_numbers( std::string s )
392 {
393 
394  int n = -1;
395  int fpos = 0;
396  double d = 0;
397  MCNPError result = MCNP_SUCCESS;
398 
399  while( result != DONE )
400  {
401  result = next_number( s, d, fpos );
402  n++;
403  }
404 
405  return n;
406 }
407 
408 MCNPError read_numbers( std::string s, int n, std::vector< double >& x )
409 {
410 
411  MCNPError result;
412  int fpos = 0;
413  double d;
414 
415  for( int i = 0; i < n; i++ )
416  {
417  result = next_number( s, d, fpos );
418  if( result == MCNP_FAILURE ) return MCNP_FAILURE;
419  x.push_back( d );
420  }
421 
422  return MCNP_SUCCESS;
423 }
424 
425 MCNPError McnpData::transform_point( double* p, double* r, int csys, double* rmat )
426 {
427 
428  double q[3];
429 
430  // Apply the rotation matrix
431  for( unsigned int i = 0; i < 3; i++ )
432  {
433  q[i] = p[0] * rmat[4 * i] + p[1] * rmat[4 * i + 1] + p[2] * rmat[4 * i + 2] + rmat[4 * i + 3];
434  }
435 
436  // Transform coordinate system
437  switch( csys )
438  {
439  case CARTESIAN:
440  r[0] = q[0];
441  r[1] = q[1];
442  r[2] = q[2]; // x, y, z
443  break;
444  case CYLINDRICAL:
445  r[0] = sqrt( q[0] * q[0] + q[1] * q[1] ); // r
446  r[1] = q[2]; // z
447  r[2] = c2pi * ( atan2( q[1], q[0] ) ); // theta (in rotations)
448  break;
449  case SPHERICAL:
450  return MCNP_FAILURE;
451  // break;
452  default:
453  return MCNP_FAILURE;
454  // break;
455  }
456 
457  return MCNP_SUCCESS;
458 }