 |
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 <cstdlib>
4 #include "mcnpmit.hpp"
5 #include "moab/CartVect.hpp"
6 #include <cmath>
7
8 moab::Interface* mb_instance();
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
20 McnpData::McnpData()
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
40 McnpData::~McnpData()
41 {
42
43 // Vertices and elements
44 MCNP_vertices.clear();
45 }
46
47 // Setting and retrieving coordinate sysem
48 MCNPError McnpData::set_coord_system( int k )
49 {
50 coord_system = k;
51 return MCNP_SUCCESS;
52 }
53 int McnpData::get_coord_system()
54 {
55 return coord_system;
56 }
57
58 // Setting and retrieving roation matrix
59 MCNPError McnpData::set_rotation_matrix( double r[16] )
60 {
61 for( int i = 0; i < 16; i++ )
62 {
63 rotation_matrix[i] = r[i];
64 }
65 return MCNP_SUCCESS;
66 }
67 double* McnpData::get_rotation_matrix()
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 }
78 std::string McnpData::get_filename()
79 {
80 return MCNP_filename;
81 }
82
83 // Reading the MCNP file
84 MCNPError McnpData::read_mcnpfile( bool skip_mesh )
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();
238 vert_handles.clear();
239 MCNP_elems.clear();
240 return MCNP_SUCCESS;
241 }
242
243 MCNPError McnpData::read_coord_system( std::string s )
244 {
245
246 if( ( s.find( "Box" ) < 100 ) || ( s.find( "xyz" ) < 100 ) )
247 coord_system = CARTESIAN;
248 else if( s.find( "Cyl" ) < 100 )
249 coord_system = CYLINDRICAL;
250 else if( s.find( "Sph" ) < 100 )
251 coord_system = SPHERICAL;
252 else
253 return MCNP_FAILURE;
254
255 return MCNP_SUCCESS;
256 }
257
258 MCNPError McnpData::read_rotation_matrix( std::string s, int i )
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
326 MCNPError McnpData::initialize_tags()
327 {
328
329 MBI->tag_get_handle( TALLY_TAG, 1, moab::MB_TYPE_DOUBLE, tally_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
330 MBI->tag_get_handle( ERROR_TAG, 1, moab::MB_TYPE_DOUBLE, relerr_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
331
332 return MCNP_SUCCESS;
333 }
334
335 MCNPError McnpData::extract_tally_data( std::string s, moab::EntityHandle handle )
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 }
1.9.1.