 |
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 /**
2 * \class ReadCGNS
3 * \brief Template for writing a new reader in MOAB
4 *
5 */
6
7 #include "ReadCGNS.hpp"
8 #include "Internals.hpp"
9 #include "moab/Interface.hpp"
10 #include "moab/ReadUtilIface.hpp"
11 #include "moab/Range.hpp"
12 #include "moab/FileOptions.hpp"
13 #include "MBTagConventions.hpp"
14 #include "MBParallelConventions.h"
15 #include "moab/CN.hpp"
16
17 #include <cstdio>
18 #include <cassert>
19 #include <cerrno>
20 #include <map>
21 #include <set>
22
23 #include <iostream>
24 #include <cmath>
25
26 namespace moab
27 {
28
29 ReaderIface* ReadCGNS::factory( Interface* iface )
30 {
31 return new ReadCGNS( iface );
32 }
33
34 ReadCGNS::ReadCGNS( Interface* impl ) : fileName( NULL ), mesh_dim( 0 ), mbImpl( impl ), globalId( 0 ), boundary( 0 )
35 {
36 mbImpl->query_interface( readMeshIface );
37 }
38
39 ReadCGNS::~ReadCGNS()
40 {
41 if( readMeshIface )
42 {
43 mbImpl->release_interface( readMeshIface );
44 readMeshIface = 0;
45 }
46 }
47
48 ErrorCode ReadCGNS::read_tag_values( const char* /* file_name */,
49 const char* /* tag_name */,
50 const FileOptions& /* opts */,
51 std::vector< int >& /* tag_values_out */,
52 const SubsetList* /* subset_list */ )
53 {
54 return MB_NOT_IMPLEMENTED;
55 }
56
57 ErrorCode ReadCGNS::load_file( const char* filename,
58 const EntityHandle* /*file_set*/,
59 const FileOptions& opts,
60 const ReaderIface::SubsetList* subset_list,
61 const Tag* file_id_tag )
62 {
63 int num_material_sets = 0;
64 const int* material_set_list = 0;
65
66 if( subset_list )
67 {
68 if( subset_list->tag_list_length > 1 && !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
69 {
70 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "CGNS supports subset read only by material ID" );
71 }
72 material_set_list = subset_list->tag_list[0].tag_values;
73 num_material_sets = subset_list->tag_list[0].num_tag_values;
74 }
75
76 ErrorCode result;
77
78 geomSets.clear();
79 globalId = mbImpl->globalId_tag();
80
81 // Create set for more convenient check for material set ids
82 std::set< int > blocks;
83 for( const int* mat_set_end = material_set_list + num_material_sets; material_set_list != mat_set_end;
84 ++material_set_list )
85 blocks.insert( *material_set_list );
86
87 // Map of ID->handle for nodes
88 std::map< long, EntityHandle > node_id_map;
89
90 // Save filename to member variable so we don't need to pass as an argument
91 // to called functions
92 fileName = filename;
93
94 // Process options; see src/FileOptions.hpp for API for FileOptions class, and
95 // doc/metadata_info.doc for a description of various options used by some of the readers in
96 // MOAB
97 result = process_options( opts );MB_CHK_SET_ERR( result, fileName << ": problem reading options" );
98
99 // Open file
100 int filePtr = 0;
101
102 cg_open( filename, CG_MODE_READ, &filePtr );
103
104 if( filePtr <= 0 )
105 {
106 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, fileName << ": fopen returned error" );
107 }
108
109 // Read number of verts, elements, sets
110 long num_verts = 0, num_elems = 0, num_sets = 0;
111 int num_bases = 0, num_zones = 0, num_sections = 0;
112
113 char zoneName[128];
114 cgsize_t size[3];
115
116 mesh_dim = 3; // Default to 3D
117
118 // Read number of bases;
119 cg_nbases( filePtr, &num_bases );
120
121 if( num_bases > 1 )
122 {
123 MB_SET_ERR( MB_NOT_IMPLEMENTED, fileName << ": support for number of bases > 1 not implemented" );
124 }
125
126 for( int indexBase = 1; indexBase <= num_bases; ++indexBase )
127 {
128 // Get the number of zones/blocks in current base.
129 cg_nzones( filePtr, indexBase, &num_zones );
130
131 if( num_zones > 1 )
132 {
133 MB_SET_ERR( MB_NOT_IMPLEMENTED, fileName << ": support for number of zones > 1 not implemented" );
134 }
135
136 for( int indexZone = 1; indexZone <= num_zones; ++indexZone )
137 {
138 // Get zone name and size.
139 cg_zone_read( filePtr, indexBase, indexZone, zoneName, size );
140
141 // Read number of sections/Parts in current zone.
142 cg_nsections( filePtr, indexBase, indexZone, &num_sections );
143
144 num_verts = size[0];
145 num_elems = size[1];
146 num_sets = num_sections;
147
148 std::cout << "\nnumber of nodes = " << num_verts;
149 std::cout << "\nnumber of elems = " << num_elems;
150 std::cout << "\nnumber of parts = " << num_sets << std::endl;
151
152 // //////////////////////////////////
153 // Read Nodes
154
155 // Allocate nodes; these are allocated in one shot, get contiguous handles starting with
156 // start_handle, and the reader is passed back double*'s pointing to MOAB's native
157 // storage for vertex coordinates for those verts
158 std::vector< double* > coord_arrays;
159 EntityHandle handle = 0;
160 result = readMeshIface->get_node_coords( 3, num_verts, MB_START_ID, handle, coord_arrays );MB_CHK_SET_ERR( result, fileName << ": Trouble reading vertices" );
161
162 // Fill in vertex coordinate arrays
163 cgsize_t beginPos = 1, endPos = num_verts;
164
165 // Read nodes coordinates.
166 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateX", RealDouble, &beginPos, &endPos,
167 coord_arrays[0] );
168 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateY", RealDouble, &beginPos, &endPos,
169 coord_arrays[1] );
170 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateZ", RealDouble, &beginPos, &endPos,
171 coord_arrays[2] );
172
173 // CGNS seems to always include the Z component, even if the mesh is 2D.
174 // Check if Z is zero and determine mesh dimension.
175 // Also create the node_id_map data.
176 double sumZcoord = 0.0;
177 double eps = 1.0e-12;
178 for( long i = 0; i < num_verts; ++i, ++handle )
179 {
180 int index = i + 1;
181
182 node_id_map.insert( std::pair< long, EntityHandle >( index, handle ) ).second;
183
184 sumZcoord += *( coord_arrays[2] + i );
185 }
186 if( std::abs( sumZcoord ) <= eps ) mesh_dim = 2;
187
188 // Create reverse map from handle to id
189 std::vector< int > ids( num_verts );
190 std::vector< int >::iterator id_iter = ids.begin();
191 std::vector< EntityHandle > handles( num_verts );
192 std::vector< EntityHandle >::iterator h_iter = handles.begin();
193 for( std::map< long, EntityHandle >::iterator i = node_id_map.begin(); i != node_id_map.end();
194 ++i, ++id_iter, ++h_iter )
195 {
196 *id_iter = i->first;
197 *h_iter = i->second;
198 }
199 // Store IDs in tags
200 result = mbImpl->tag_set_data( globalId, &handles[0], num_verts, &ids[0] );
201 if( MB_SUCCESS != result ) return result;
202 if( file_id_tag )
203 {
204 result = mbImpl->tag_set_data( *file_id_tag, &handles[0], num_verts, &ids[0] );
205 if( MB_SUCCESS != result ) return result;
206 }
207 ids.clear();
208 handles.clear();
209
210 // //////////////////////////////////
211 // Read elements data
212
213 EntityType ent_type;
214
215 long section_offset = 0;
216
217 // Define which mesh parts are volume families.
218 // Mesh parts with volumeID[X] = 0 are boundary parts.
219 std::vector< int > volumeID( num_sections, 0 );
220
221 for( int section = 0; section < num_sections; ++section )
222 {
223 ElementType_t elemsType;
224 int iparent_flag, nbndry;
225 char sectionName[128];
226 int verts_per_elem;
227
228 int cgSection = section + 1;
229
230 cg_section_read( filePtr, indexBase, indexZone, cgSection, sectionName, &elemsType, &beginPos, &endPos,
231 &nbndry, &iparent_flag );
232
233 size_t section_size = endPos - beginPos + 1;
234
235 // Read element description in current section
236
237 switch( elemsType )
238 {
239 case BAR_2:
240 ent_type = MBEDGE;
241 verts_per_elem = 2;
242 break;
243 case TRI_3:
244 ent_type = MBTRI;
245 verts_per_elem = 3;
246 if( mesh_dim == 2 ) volumeID[section] = 1;
247 break;
248 case QUAD_4:
249 ent_type = MBQUAD;
250 verts_per_elem = 4;
251 if( mesh_dim == 2 ) volumeID[section] = 1;
252 break;
253 case TETRA_4:
254 ent_type = MBTET;
255 verts_per_elem = 4;
256 if( mesh_dim == 3 ) volumeID[section] = 1;
257 break;
258 case PYRA_5:
259 ent_type = MBPYRAMID;
260 verts_per_elem = 5;
261 if( mesh_dim == 3 ) volumeID[section] = 1;
262 break;
263 case PENTA_6:
264 ent_type = MBPRISM;
265 verts_per_elem = 6;
266 if( mesh_dim == 3 ) volumeID[section] = 1;
267 break;
268 case HEXA_8:
269 ent_type = MBHEX;
270 verts_per_elem = 8;
271 if( mesh_dim == 3 ) volumeID[section] = 1;
272 break;
273 case MIXED:
274 ent_type = MBMAXTYPE;
275 verts_per_elem = 0;
276 break;
277 default:
278 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" );
279 }
280
281 if( elemsType == TETRA_4 || elemsType == PYRA_5 || elemsType == PENTA_6 || elemsType == HEXA_8 ||
282 elemsType == TRI_3 || elemsType == QUAD_4 || ( ( elemsType == BAR_2 ) && mesh_dim == 2 ) )
283 {
284 // Read connectivity into conn_array directly
285
286 cgsize_t iparentdata;
287 cgsize_t connDataSize;
288
289 // Get number of entries on the connectivity list for this section
290 cg_ElementDataSize( filePtr, indexBase, indexZone, cgSection, &connDataSize );
291
292 // Need a temporary vector to later cast to conn_array.
293 std::vector< cgsize_t > elemNodes( connDataSize );
294
295 cg_elements_read( filePtr, indexBase, indexZone, cgSection, &elemNodes[0], &iparentdata );
296
297 // //////////////////////////////////
298 // Create elements, sets and tags
299
300 create_elements( sectionName, file_id_tag, ent_type, verts_per_elem, section_offset, section_size,
301 elemNodes );
302 } // Homogeneous mesh type
303 else if( elemsType == MIXED )
304 {
305 // We must first sort all elements connectivities into continuous vectors
306
307 cgsize_t connDataSize;
308 cgsize_t iparentdata;
309
310 cg_ElementDataSize( filePtr, indexBase, indexZone, cgSection, &connDataSize );
311
312 std::vector< cgsize_t > elemNodes( connDataSize );
313
314 cg_elements_read( filePtr, indexBase, indexZone, cgSection, &elemNodes[0], &iparentdata );
315
316 std::vector< cgsize_t > elemsConn_EDGE;
317 std::vector< cgsize_t > elemsConn_TRI, elemsConn_QUAD;
318 std::vector< cgsize_t > elemsConn_TET, elemsConn_PYRA, elemsConn_PRISM, elemsConn_HEX;
319 cgsize_t count_EDGE, count_TRI, count_QUAD;
320 cgsize_t count_TET, count_PYRA, count_PRISM, count_HEX;
321
322 // First, get elements count for current section
323
324 count_EDGE = count_TRI = count_QUAD = 0;
325 count_TET = count_PYRA = count_PRISM = count_HEX = 0;
326
327 int connIndex = 0;
328 for( int i = beginPos; i <= endPos; i++ )
329 {
330 elemsType = ElementType_t( elemNodes[connIndex] );
331
332 // Get current cell node count.
333 cg_npe( elemsType, &verts_per_elem );
334
335 switch( elemsType )
336 {
337 case BAR_2:
338 count_EDGE += 1;
339 break;
340 case TRI_3:
341 count_TRI += 1;
342 break;
343 case QUAD_4:
344 count_QUAD += 1;
345 break;
346 case TETRA_4:
347 count_TET += 1;
348 break;
349 case PYRA_5:
350 count_PYRA += 1;
351 break;
352 case PENTA_6:
353 count_PRISM += 1;
354 break;
355 case HEXA_8:
356 count_HEX += 1;
357 break;
358 default:
359 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" );
360 }
361
362 connIndex += ( verts_per_elem + 1 ); // Add one to skip next element descriptor
363 }
364
365 if( count_EDGE > 0 ) elemsConn_EDGE.resize( count_EDGE * 2 );
366 if( count_TRI > 0 ) elemsConn_TRI.resize( count_TRI * 3 );
367 if( count_QUAD > 0 ) elemsConn_QUAD.resize( count_QUAD * 4 );
368 if( count_TET > 0 ) elemsConn_TET.resize( count_TET * 4 );
369 if( count_PYRA > 0 ) elemsConn_PYRA.resize( count_PYRA * 5 );
370 if( count_PRISM > 0 ) elemsConn_PRISM.resize( count_PRISM * 6 );
371 if( count_HEX > 0 ) elemsConn_HEX.resize( count_HEX * 8 );
372
373 // Grab mixed section elements connectivity
374
375 int idx_edge, idx_tri, idx_quad;
376 int idx_tet, idx_pyra, idx_prism, idx_hex;
377 idx_edge = idx_tri = idx_quad = 0;
378 idx_tet = idx_pyra = idx_prism = idx_hex = 0;
379
380 connIndex = 0;
381 for( int i = beginPos; i <= endPos; i++ )
382 {
383 elemsType = ElementType_t( elemNodes[connIndex] );
384
385 // Get current cell node count.
386 cg_npe( elemsType, &verts_per_elem );
387
388 switch( elemsType )
389 {
390 case BAR_2:
391 for( int j = 0; j < 2; ++j )
392 elemsConn_EDGE[idx_edge + j] = elemNodes[connIndex + j + 1];
393 idx_edge += 2;
394 break;
395 case TRI_3:
396 for( int j = 0; j < 3; ++j )
397 elemsConn_TRI[idx_tri + j] = elemNodes[connIndex + j + 1];
398 idx_tri += 3;
399 break;
400 case QUAD_4:
401 for( int j = 0; j < 4; ++j )
402 elemsConn_QUAD[idx_quad + j] = elemNodes[connIndex + j + 1];
403 idx_quad += 4;
404 break;
405 case TETRA_4:
406 for( int j = 0; j < 4; ++j )
407 elemsConn_TET[idx_tet + j] = elemNodes[connIndex + j + 1];
408 idx_tet += 4;
409 break;
410 case PYRA_5:
411 for( int j = 0; j < 5; ++j )
412 elemsConn_PYRA[idx_pyra + j] = elemNodes[connIndex + j + 1];
413 idx_pyra += 5;
414 break;
415 case PENTA_6:
416 for( int j = 0; j < 6; ++j )
417 elemsConn_PRISM[idx_prism + j] = elemNodes[connIndex + j + 1];
418 idx_prism += 6;
419 break;
420 case HEXA_8:
421 for( int j = 0; j < 8; ++j )
422 elemsConn_HEX[idx_hex + j] = elemNodes[connIndex + j + 1];
423 idx_hex += 8;
424 break;
425 default:
426 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" );
427 }
428
429 connIndex += ( verts_per_elem + 1 ); // Add one to skip next element descriptor
430 }
431
432 // //////////////////////////////////
433 // Create elements, sets and tags
434
435 if( count_EDGE > 0 )
436 create_elements( sectionName, file_id_tag, MBEDGE, 2, section_offset, count_EDGE,
437 elemsConn_EDGE );
438
439 if( count_TRI > 0 )
440 create_elements( sectionName, file_id_tag, MBTRI, 3, section_offset, count_TRI, elemsConn_TRI );
441
442 if( count_QUAD > 0 )
443 create_elements( sectionName, file_id_tag, MBQUAD, 4, section_offset, count_QUAD,
444 elemsConn_QUAD );
445
446 if( count_TET > 0 )
447 create_elements( sectionName, file_id_tag, MBTET, 4, section_offset, count_TET, elemsConn_TET );
448
449 if( count_PYRA > 0 )
450 create_elements( sectionName, file_id_tag, MBPYRAMID, 5, section_offset, count_PYRA,
451 elemsConn_PYRA );
452
453 if( count_PRISM > 0 )
454 create_elements( sectionName, file_id_tag, MBPRISM, 6, section_offset, count_PRISM,
455 elemsConn_PRISM );
456
457 if( count_HEX > 0 )
458 create_elements( sectionName, file_id_tag, MBHEX, 8, section_offset, count_HEX, elemsConn_HEX );
459 } // Mixed mesh type
460 } // num_sections
461
462 cg_close( filePtr );
463
464 return result;
465 } // indexZone for
466 } // indexBase for
467
468 return MB_SUCCESS;
469 }
470
471 ErrorCode ReadCGNS::create_elements( char* sectionName,
472 const Tag* file_id_tag,
473 const EntityType& ent_type,
474 const int& verts_per_elem,
475 long& section_offset,
476 int elems_count,
477 const std::vector< cgsize_t >& elemsConn )
478 {
479 ErrorCode result;
480
481 // Create the element sequence; passes back a pointer to the internal storage for connectivity
482 // and the starting entity handle
483 EntityHandle* conn_array;
484 EntityHandle handle = 0;
485
486 result = readMeshIface->get_element_connect( elems_count, verts_per_elem, ent_type, 1, handle, conn_array );MB_CHK_SET_ERR( result, fileName << ": Trouble reading elements" );
487
488 if( sizeof( EntityHandle ) == sizeof( cgsize_t ) )
489 {
490 memcpy( conn_array, &elemsConn[0], elemsConn.size() * sizeof( EntityHandle ) );
491 }
492 else
493 { // if CGNS is compiled without 64bit enabled
494 std::vector< EntityHandle > elemsConnTwin( elemsConn.size(), 0 );
495 for( int i = 0; i < elemsConn.size(); i++ )
496 {
497 elemsConnTwin[i] = static_cast< EntityHandle >( elemsConn[i] );
498 }
499 memcpy( conn_array, &elemsConnTwin[0], elemsConnTwin.size() * sizeof( EntityHandle ) );
500 }
501
502 // Notify MOAB of the new elements
503 result = readMeshIface->update_adjacencies( handle, elems_count, verts_per_elem, conn_array );
504 if( MB_SUCCESS != result ) return result;
505
506 // //////////////////////////////////
507 // Create sets and tags
508
509 Range elements( handle, handle + elems_count - 1 );
510
511 // Store element IDs
512
513 std::vector< int > id_list( elems_count );
514
515 // Add 1 to offset id to 1-based numbering
516 for( cgsize_t i = 0; i < elems_count; ++i )
517 id_list[i] = i + 1 + section_offset;
518 section_offset += elems_count;
519
520 create_sets( sectionName, file_id_tag, ent_type, elements, id_list, 0 );
521
522 return MB_SUCCESS;
523 }
524
525 ErrorCode ReadCGNS::create_sets( char* sectionName,
526 const Tag* file_id_tag,
527 EntityType /*element_type*/,
528 const Range& elements,
529 const std::vector< int >& set_ids,
530 int /*set_type*/ )
531 {
532 ErrorCode result;
533
534 result = mbImpl->tag_set_data( globalId, elements, &set_ids[0] );
535 if( MB_SUCCESS != result ) return result;
536
537 if( file_id_tag )
538 {
539 result = mbImpl->tag_set_data( *file_id_tag, elements, &set_ids[0] );
540 if( MB_SUCCESS != result ) return result;
541 }
542
543 EntityHandle set_handle;
544
545 Tag tag_handle;
546
547 const char* setName = sectionName;
548
549 mbImpl->tag_get_handle( setName, 1, MB_TYPE_INTEGER, tag_handle, MB_TAG_SPARSE | MB_TAG_CREAT );
550
551 // Create set
552 result = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( result, fileName << ": Trouble creating set" );
553
554 //// Add dummy values to current set
555 // std::vector<int> tags(set_ids.size(), 1);
556 // result = mbImpl->tag_set_data(tag_handle, elements, &tags[0]);
557 // if (MB_SUCCESS != result) return result;
558
559 // Add them to the set
560 result = mbImpl->add_entities( set_handle, elements );MB_CHK_SET_ERR( result, fileName << ": Trouble putting entities in set" );
561
562 return MB_SUCCESS;
563 }
564
565 ErrorCode ReadCGNS::process_options( const FileOptions& opts )
566 {
567 // Mark all options seen, to avoid compile warning on unused variable
568 opts.mark_all_seen();
569
570 return MB_SUCCESS;
571 }
572
573 } // namespace moab
1.9.1.