MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * \class ReadCGNS 00003 * \brief Template for writing a new reader in MOAB 00004 * 00005 */ 00006 00007 #include "ReadCGNS.hpp" 00008 #include "Internals.hpp" 00009 #include "moab/Interface.hpp" 00010 #include "moab/ReadUtilIface.hpp" 00011 #include "moab/Range.hpp" 00012 #include "moab/FileOptions.hpp" 00013 #include "MBTagConventions.hpp" 00014 #include "MBParallelConventions.h" 00015 #include "moab/CN.hpp" 00016 00017 #include <cstdio> 00018 #include <assert.h> 00019 #include <errno.h> 00020 #include <map> 00021 #include <set> 00022 00023 #include <iostream> 00024 #include <cmath> 00025 00026 namespace moab 00027 { 00028 00029 ReaderIface* ReadCGNS::factory( Interface* iface ) 00030 { 00031 return new ReadCGNS( iface ); 00032 } 00033 00034 ReadCGNS::ReadCGNS( Interface* impl ) : fileName( NULL ), mesh_dim( 0 ), mbImpl( impl ), globalId( 0 ), boundary( 0 ) 00035 { 00036 mbImpl->query_interface( readMeshIface ); 00037 } 00038 00039 ReadCGNS::~ReadCGNS() 00040 { 00041 if( readMeshIface ) 00042 { 00043 mbImpl->release_interface( readMeshIface ); 00044 readMeshIface = 0; 00045 } 00046 } 00047 00048 ErrorCode ReadCGNS::read_tag_values( const char* /* file_name */, const char* /* tag_name */, 00049 const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */, 00050 const SubsetList* /* subset_list */ ) 00051 { 00052 return MB_NOT_IMPLEMENTED; 00053 } 00054 00055 ErrorCode ReadCGNS::load_file( const char* filename, const EntityHandle* /*file_set*/, const FileOptions& opts, 00056 const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag ) 00057 { 00058 int num_material_sets = 0; 00059 const int* material_set_list = 0; 00060 00061 if( subset_list ) 00062 { 00063 if( subset_list->tag_list_length > 1 && !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) ) 00064 { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "CGNS supports subset read only by material ID" ); } 00065 material_set_list = subset_list->tag_list[0].tag_values; 00066 num_material_sets = subset_list->tag_list[0].num_tag_values; 00067 } 00068 00069 ErrorCode result; 00070 00071 geomSets.clear(); 00072 globalId = mbImpl->globalId_tag(); 00073 00074 // Create set for more convenient check for material set ids 00075 std::set< int > blocks; 00076 for( const int* mat_set_end = material_set_list + num_material_sets; material_set_list != mat_set_end; 00077 ++material_set_list ) 00078 blocks.insert( *material_set_list ); 00079 00080 // Map of ID->handle for nodes 00081 std::map< long, EntityHandle > node_id_map; 00082 00083 // Save filename to member variable so we don't need to pass as an argument 00084 // to called functions 00085 fileName = filename; 00086 00087 // Process options; see src/FileOptions.hpp for API for FileOptions class, and 00088 // doc/metadata_info.doc for a description of various options used by some of the readers in 00089 // MOAB 00090 result = process_options( opts );MB_CHK_SET_ERR( result, fileName << ": problem reading options" ); 00091 00092 // Open file 00093 int filePtr = 0; 00094 00095 cg_open( filename, CG_MODE_READ, &filePtr ); 00096 00097 if( filePtr <= 0 ) { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, fileName << ": fopen returned error" ); } 00098 00099 // Read number of verts, elements, sets 00100 long num_verts = 0, num_elems = 0, num_sets = 0; 00101 int num_bases = 0, num_zones = 0, num_sections = 0; 00102 00103 char zoneName[128]; 00104 cgsize_t size[3]; 00105 00106 mesh_dim = 3; // Default to 3D 00107 00108 // Read number of bases; 00109 cg_nbases( filePtr, &num_bases ); 00110 00111 if( num_bases > 1 ) 00112 { MB_SET_ERR( MB_NOT_IMPLEMENTED, fileName << ": support for number of bases > 1 not implemented" ); } 00113 00114 for( int indexBase = 1; indexBase <= num_bases; ++indexBase ) 00115 { 00116 // Get the number of zones/blocks in current base. 00117 cg_nzones( filePtr, indexBase, &num_zones ); 00118 00119 if( num_zones > 1 ) 00120 { MB_SET_ERR( MB_NOT_IMPLEMENTED, fileName << ": support for number of zones > 1 not implemented" ); } 00121 00122 for( int indexZone = 1; indexZone <= num_zones; ++indexZone ) 00123 { 00124 // Get zone name and size. 00125 cg_zone_read( filePtr, indexBase, indexZone, zoneName, size ); 00126 00127 // Read number of sections/Parts in current zone. 00128 cg_nsections( filePtr, indexBase, indexZone, &num_sections ); 00129 00130 num_verts = size[0]; 00131 num_elems = size[1]; 00132 num_sets = num_sections; 00133 00134 std::cout << "\nnumber of nodes = " << num_verts; 00135 std::cout << "\nnumber of elems = " << num_elems; 00136 std::cout << "\nnumber of parts = " << num_sets << std::endl; 00137 00138 // ////////////////////////////////// 00139 // Read Nodes 00140 00141 // Allocate nodes; these are allocated in one shot, get contiguous handles starting with 00142 // start_handle, and the reader is passed back double*'s pointing to MOAB's native 00143 // storage for vertex coordinates for those verts 00144 std::vector< double* > coord_arrays; 00145 EntityHandle handle = 0; 00146 result = readMeshIface->get_node_coords( 3, num_verts, MB_START_ID, handle, coord_arrays );MB_CHK_SET_ERR( result, fileName << ": Trouble reading vertices" ); 00147 00148 // Fill in vertex coordinate arrays 00149 cgsize_t beginPos = 1, endPos = num_verts; 00150 00151 // Read nodes coordinates. 00152 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateX", RealDouble, &beginPos, &endPos, 00153 coord_arrays[0] ); 00154 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateY", RealDouble, &beginPos, &endPos, 00155 coord_arrays[1] ); 00156 cg_coord_read( filePtr, indexBase, indexZone, "CoordinateZ", RealDouble, &beginPos, &endPos, 00157 coord_arrays[2] ); 00158 00159 // CGNS seems to always include the Z component, even if the mesh is 2D. 00160 // Check if Z is zero and determine mesh dimension. 00161 // Also create the node_id_map data. 00162 double sumZcoord = 0.0; 00163 double eps = 1.0e-12; 00164 for( long i = 0; i < num_verts; ++i, ++handle ) 00165 { 00166 int index = i + 1; 00167 00168 node_id_map.insert( std::pair< long, EntityHandle >( index, handle ) ).second; 00169 00170 sumZcoord += *( coord_arrays[2] + i ); 00171 } 00172 if( std::abs( sumZcoord ) <= eps ) mesh_dim = 2; 00173 00174 // Create reverse map from handle to id 00175 std::vector< int > ids( num_verts ); 00176 std::vector< int >::iterator id_iter = ids.begin(); 00177 std::vector< EntityHandle > handles( num_verts ); 00178 std::vector< EntityHandle >::iterator h_iter = handles.begin(); 00179 for( std::map< long, EntityHandle >::iterator i = node_id_map.begin(); i != node_id_map.end(); 00180 ++i, ++id_iter, ++h_iter ) 00181 { 00182 *id_iter = i->first; 00183 *h_iter = i->second; 00184 } 00185 // Store IDs in tags 00186 result = mbImpl->tag_set_data( globalId, &handles[0], num_verts, &ids[0] ); 00187 if( MB_SUCCESS != result ) return result; 00188 if( file_id_tag ) 00189 { 00190 result = mbImpl->tag_set_data( *file_id_tag, &handles[0], num_verts, &ids[0] ); 00191 if( MB_SUCCESS != result ) return result; 00192 } 00193 ids.clear(); 00194 handles.clear(); 00195 00196 // ////////////////////////////////// 00197 // Read elements data 00198 00199 EntityType ent_type; 00200 00201 long section_offset = 0; 00202 00203 // Define which mesh parts are volume families. 00204 // Mesh parts with volumeID[X] = 0 are boundary parts. 00205 std::vector< int > volumeID( num_sections, 0 ); 00206 00207 for( int section = 0; section < num_sections; ++section ) 00208 { 00209 ElementType_t elemsType; 00210 int iparent_flag, nbndry; 00211 char sectionName[128]; 00212 int verts_per_elem; 00213 00214 int cgSection = section + 1; 00215 00216 cg_section_read( filePtr, indexBase, indexZone, cgSection, sectionName, &elemsType, &beginPos, &endPos, 00217 &nbndry, &iparent_flag ); 00218 00219 size_t section_size = endPos - beginPos + 1; 00220 00221 // Read element description in current section 00222 00223 switch( elemsType ) 00224 { 00225 case BAR_2: 00226 ent_type = MBEDGE; 00227 verts_per_elem = 2; 00228 break; 00229 case TRI_3: 00230 ent_type = MBTRI; 00231 verts_per_elem = 3; 00232 if( mesh_dim == 2 ) volumeID[section] = 1; 00233 break; 00234 case QUAD_4: 00235 ent_type = MBQUAD; 00236 verts_per_elem = 4; 00237 if( mesh_dim == 2 ) volumeID[section] = 1; 00238 break; 00239 case TETRA_4: 00240 ent_type = MBTET; 00241 verts_per_elem = 4; 00242 if( mesh_dim == 3 ) volumeID[section] = 1; 00243 break; 00244 case PYRA_5: 00245 ent_type = MBPYRAMID; 00246 verts_per_elem = 5; 00247 if( mesh_dim == 3 ) volumeID[section] = 1; 00248 break; 00249 case PENTA_6: 00250 ent_type = MBPRISM; 00251 verts_per_elem = 6; 00252 if( mesh_dim == 3 ) volumeID[section] = 1; 00253 break; 00254 case HEXA_8: 00255 ent_type = MBHEX; 00256 verts_per_elem = 8; 00257 if( mesh_dim == 3 ) volumeID[section] = 1; 00258 break; 00259 case MIXED: 00260 ent_type = MBMAXTYPE; 00261 verts_per_elem = 0; 00262 break; 00263 default: 00264 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" ); 00265 } 00266 00267 if( elemsType == TETRA_4 || elemsType == PYRA_5 || elemsType == PENTA_6 || elemsType == HEXA_8 || 00268 elemsType == TRI_3 || elemsType == QUAD_4 || ( ( elemsType == BAR_2 ) && mesh_dim == 2 ) ) 00269 { 00270 // Read connectivity into conn_array directly 00271 00272 cgsize_t iparentdata; 00273 cgsize_t connDataSize; 00274 00275 // Get number of entries on the connectivity list for this section 00276 cg_ElementDataSize( filePtr, indexBase, indexZone, cgSection, &connDataSize ); 00277 00278 // Need a temporary vector to later cast to conn_array. 00279 std::vector< cgsize_t > elemNodes( connDataSize ); 00280 00281 cg_elements_read( filePtr, indexBase, indexZone, cgSection, &elemNodes[0], &iparentdata ); 00282 00283 // ////////////////////////////////// 00284 // Create elements, sets and tags 00285 00286 create_elements( sectionName, file_id_tag, ent_type, verts_per_elem, section_offset, section_size, 00287 elemNodes ); 00288 } // Homogeneous mesh type 00289 else if( elemsType == MIXED ) 00290 { 00291 // We must first sort all elements connectivities into continuous vectors 00292 00293 cgsize_t connDataSize; 00294 cgsize_t iparentdata; 00295 00296 cg_ElementDataSize( filePtr, indexBase, indexZone, cgSection, &connDataSize ); 00297 00298 std::vector< cgsize_t > elemNodes( connDataSize ); 00299 00300 cg_elements_read( filePtr, indexBase, indexZone, cgSection, &elemNodes[0], &iparentdata ); 00301 00302 std::vector< cgsize_t > elemsConn_EDGE; 00303 std::vector< cgsize_t > elemsConn_TRI, elemsConn_QUAD; 00304 std::vector< cgsize_t > elemsConn_TET, elemsConn_PYRA, elemsConn_PRISM, elemsConn_HEX; 00305 cgsize_t count_EDGE, count_TRI, count_QUAD; 00306 cgsize_t count_TET, count_PYRA, count_PRISM, count_HEX; 00307 00308 // First, get elements count for current section 00309 00310 count_EDGE = count_TRI = count_QUAD = 0; 00311 count_TET = count_PYRA = count_PRISM = count_HEX = 0; 00312 00313 int connIndex = 0; 00314 for( int i = beginPos; i <= endPos; i++ ) 00315 { 00316 elemsType = ElementType_t( elemNodes[connIndex] ); 00317 00318 // Get current cell node count. 00319 cg_npe( elemsType, &verts_per_elem ); 00320 00321 switch( elemsType ) 00322 { 00323 case BAR_2: 00324 count_EDGE += 1; 00325 break; 00326 case TRI_3: 00327 count_TRI += 1; 00328 break; 00329 case QUAD_4: 00330 count_QUAD += 1; 00331 break; 00332 case TETRA_4: 00333 count_TET += 1; 00334 break; 00335 case PYRA_5: 00336 count_PYRA += 1; 00337 break; 00338 case PENTA_6: 00339 count_PRISM += 1; 00340 break; 00341 case HEXA_8: 00342 count_HEX += 1; 00343 break; 00344 default: 00345 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" ); 00346 } 00347 00348 connIndex += ( verts_per_elem + 1 ); // Add one to skip next element descriptor 00349 } 00350 00351 if( count_EDGE > 0 ) elemsConn_EDGE.resize( count_EDGE * 2 ); 00352 if( count_TRI > 0 ) elemsConn_TRI.resize( count_TRI * 3 ); 00353 if( count_QUAD > 0 ) elemsConn_QUAD.resize( count_QUAD * 4 ); 00354 if( count_TET > 0 ) elemsConn_TET.resize( count_TET * 4 ); 00355 if( count_PYRA > 0 ) elemsConn_PYRA.resize( count_PYRA * 5 ); 00356 if( count_PRISM > 0 ) elemsConn_PRISM.resize( count_PRISM * 6 ); 00357 if( count_HEX > 0 ) elemsConn_HEX.resize( count_HEX * 8 ); 00358 00359 // Grab mixed section elements connectivity 00360 00361 int idx_edge, idx_tri, idx_quad; 00362 int idx_tet, idx_pyra, idx_prism, idx_hex; 00363 idx_edge = idx_tri = idx_quad = 0; 00364 idx_tet = idx_pyra = idx_prism = idx_hex = 0; 00365 00366 connIndex = 0; 00367 for( int i = beginPos; i <= endPos; i++ ) 00368 { 00369 elemsType = ElementType_t( elemNodes[connIndex] ); 00370 00371 // Get current cell node count. 00372 cg_npe( elemsType, &verts_per_elem ); 00373 00374 switch( elemsType ) 00375 { 00376 case BAR_2: 00377 for( int j = 0; j < 2; ++j ) 00378 elemsConn_EDGE[idx_edge + j] = elemNodes[connIndex + j + 1]; 00379 idx_edge += 2; 00380 break; 00381 case TRI_3: 00382 for( int j = 0; j < 3; ++j ) 00383 elemsConn_TRI[idx_tri + j] = elemNodes[connIndex + j + 1]; 00384 idx_tri += 3; 00385 break; 00386 case QUAD_4: 00387 for( int j = 0; j < 4; ++j ) 00388 elemsConn_QUAD[idx_quad + j] = elemNodes[connIndex + j + 1]; 00389 idx_quad += 4; 00390 break; 00391 case TETRA_4: 00392 for( int j = 0; j < 4; ++j ) 00393 elemsConn_TET[idx_tet + j] = elemNodes[connIndex + j + 1]; 00394 idx_tet += 4; 00395 break; 00396 case PYRA_5: 00397 for( int j = 0; j < 5; ++j ) 00398 elemsConn_PYRA[idx_pyra + j] = elemNodes[connIndex + j + 1]; 00399 idx_pyra += 5; 00400 break; 00401 case PENTA_6: 00402 for( int j = 0; j < 6; ++j ) 00403 elemsConn_PRISM[idx_prism + j] = elemNodes[connIndex + j + 1]; 00404 idx_prism += 6; 00405 break; 00406 case HEXA_8: 00407 for( int j = 0; j < 8; ++j ) 00408 elemsConn_HEX[idx_hex + j] = elemNodes[connIndex + j + 1]; 00409 idx_hex += 8; 00410 break; 00411 default: 00412 MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, fileName << ": Trouble determining element type" ); 00413 } 00414 00415 connIndex += ( verts_per_elem + 1 ); // Add one to skip next element descriptor 00416 } 00417 00418 // ////////////////////////////////// 00419 // Create elements, sets and tags 00420 00421 if( count_EDGE > 0 ) 00422 create_elements( sectionName, file_id_tag, MBEDGE, 2, section_offset, count_EDGE, 00423 elemsConn_EDGE ); 00424 00425 if( count_TRI > 0 ) 00426 create_elements( sectionName, file_id_tag, MBTRI, 3, section_offset, count_TRI, elemsConn_TRI ); 00427 00428 if( count_QUAD > 0 ) 00429 create_elements( sectionName, file_id_tag, MBQUAD, 4, section_offset, count_QUAD, 00430 elemsConn_QUAD ); 00431 00432 if( count_TET > 0 ) 00433 create_elements( sectionName, file_id_tag, MBTET, 4, section_offset, count_TET, elemsConn_TET ); 00434 00435 if( count_PYRA > 0 ) 00436 create_elements( sectionName, file_id_tag, MBPYRAMID, 5, section_offset, count_PYRA, 00437 elemsConn_PYRA ); 00438 00439 if( count_PRISM > 0 ) 00440 create_elements( sectionName, file_id_tag, MBPRISM, 6, section_offset, count_PRISM, 00441 elemsConn_PRISM ); 00442 00443 if( count_HEX > 0 ) 00444 create_elements( sectionName, file_id_tag, MBHEX, 8, section_offset, count_HEX, elemsConn_HEX ); 00445 } // Mixed mesh type 00446 } // num_sections 00447 00448 cg_close( filePtr ); 00449 00450 return result; 00451 } // indexZone for 00452 } // indexBase for 00453 00454 return MB_SUCCESS; 00455 } 00456 00457 ErrorCode ReadCGNS::create_elements( char* sectionName, const Tag* file_id_tag, const EntityType& ent_type, 00458 const int& verts_per_elem, long& section_offset, int elems_count, 00459 const std::vector< cgsize_t >& elemsConn ) 00460 { 00461 ErrorCode result; 00462 00463 // Create the element sequence; passes back a pointer to the internal storage for connectivity 00464 // and the starting entity handle 00465 EntityHandle* conn_array; 00466 EntityHandle handle = 0; 00467 00468 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" ); 00469 00470 memcpy( conn_array, &elemsConn[0], elemsConn.size() * sizeof( EntityHandle ) ); 00471 00472 // Notify MOAB of the new elements 00473 result = readMeshIface->update_adjacencies( handle, elems_count, verts_per_elem, conn_array ); 00474 if( MB_SUCCESS != result ) return result; 00475 00476 // ////////////////////////////////// 00477 // Create sets and tags 00478 00479 Range elements( handle, handle + elems_count - 1 ); 00480 00481 // Store element IDs 00482 00483 std::vector< int > id_list( elems_count ); 00484 00485 // Add 1 to offset id to 1-based numbering 00486 for( cgsize_t i = 0; i < elems_count; ++i ) 00487 id_list[i] = i + 1 + section_offset; 00488 section_offset += elems_count; 00489 00490 create_sets( sectionName, file_id_tag, ent_type, elements, id_list, 0 ); 00491 00492 return MB_SUCCESS; 00493 } 00494 00495 ErrorCode ReadCGNS::create_sets( char* sectionName, const Tag* file_id_tag, EntityType /*element_type*/, 00496 const Range& elements, const std::vector< int >& set_ids, int /*set_type*/ ) 00497 { 00498 ErrorCode result; 00499 00500 result = mbImpl->tag_set_data( globalId, elements, &set_ids[0] ); 00501 if( MB_SUCCESS != result ) return result; 00502 00503 if( file_id_tag ) 00504 { 00505 result = mbImpl->tag_set_data( *file_id_tag, elements, &set_ids[0] ); 00506 if( MB_SUCCESS != result ) return result; 00507 } 00508 00509 EntityHandle set_handle; 00510 00511 Tag tag_handle; 00512 00513 const char* setName = sectionName; 00514 00515 mbImpl->tag_get_handle( setName, 1, MB_TYPE_INTEGER, tag_handle, MB_TAG_SPARSE | MB_TAG_CREAT ); 00516 00517 // Create set 00518 result = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( result, fileName << ": Trouble creating set" ); 00519 00520 //// Add dummy values to current set 00521 // std::vector<int> tags(set_ids.size(), 1); 00522 // result = mbImpl->tag_set_data(tag_handle, elements, &tags[0]); 00523 // if (MB_SUCCESS != result) return result; 00524 00525 // Add them to the set 00526 result = mbImpl->add_entities( set_handle, elements );MB_CHK_SET_ERR( result, fileName << ": Trouble putting entities in set" ); 00527 00528 return MB_SUCCESS; 00529 } 00530 00531 ErrorCode ReadCGNS::process_options( const FileOptions& opts ) 00532 { 00533 // Mark all options seen, to avoid compile warning on unused variable 00534 opts.mark_all_seen(); 00535 00536 return MB_SUCCESS; 00537 } 00538 00539 } // namespace moab