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