moab
|
00001 00025 #include "ReadGmsh.hpp" 00026 #include "FileTokenizer.hpp" // for file tokenizer 00027 #include "Internals.hpp" 00028 #include "moab/Interface.hpp" 00029 #include "moab/ReadUtilIface.hpp" 00030 #include "moab/Range.hpp" 00031 #include "MBTagConventions.hpp" 00032 #include "MBParallelConventions.h" 00033 #include "moab/CN.hpp" 00034 #include "GmshUtil.hpp" 00035 00036 #include <errno.h> 00037 #include <string.h> 00038 #include <map> 00039 #include <set> 00040 00041 namespace moab { 00042 00043 ReaderIface* ReadGmsh::factory( Interface* iface ) 00044 { return new ReadGmsh(iface); } 00045 00046 ReadGmsh::ReadGmsh(Interface* impl) 00047 : mdbImpl(impl) 00048 { 00049 mdbImpl->query_interface(readMeshIface); 00050 } 00051 00052 ReadGmsh::~ReadGmsh() 00053 { 00054 if (readMeshIface) { 00055 mdbImpl->release_interface(readMeshIface); 00056 readMeshIface = 0; 00057 } 00058 } 00059 00060 00061 ErrorCode ReadGmsh::read_tag_values( const char* /* file_name */, 00062 const char* /* tag_name */, 00063 const FileOptions& /* opts */, 00064 std::vector<int>& /* tag_values_out */, 00065 const SubsetList* /* subset_list */ ) 00066 { 00067 return MB_NOT_IMPLEMENTED; 00068 } 00069 00070 00071 ErrorCode ReadGmsh::load_file( const char* filename, 00072 const EntityHandle*, 00073 const FileOptions& , 00074 const ReaderIface::SubsetList* subset_list, 00075 const Tag* file_id_tag ) 00076 { 00077 int num_material_sets = 0; 00078 const int* material_set_list = 0; 00079 int zero = 0; 00080 if (subset_list) { 00081 if (subset_list->tag_list_length > 1 && 00082 !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME) ) { 00083 readMeshIface->report_error( "GMsh supports subset read only by material ID." ); 00084 return MB_UNSUPPORTED_OPERATION; 00085 } 00086 material_set_list = subset_list->tag_list[0].tag_values; 00087 num_material_sets = subset_list->tag_list[0].num_tag_values; 00088 } 00089 00090 geomSets.clear(); 00091 ErrorCode result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, 00092 globalId, MB_TAG_DENSE|MB_TAG_CREAT, &zero); 00093 if (MB_SUCCESS != result) 00094 return result; 00095 00096 // Create set for more convienient check for material set ids 00097 std::set<int> blocks; 00098 for (const int* mat_set_end = material_set_list + num_material_sets; 00099 material_set_list != mat_set_end; ++material_set_list) 00100 blocks.insert( *material_set_list ); 00101 00102 // Map of ID->handle for nodes 00103 std::map<long,EntityHandle> node_id_map; 00104 int data_size = 8; 00105 00106 // Open file and hand off pointer to tokenizer 00107 FILE* file_ptr = fopen( filename, "r" ); 00108 if (!file_ptr) 00109 { 00110 readMeshIface->report_error( "%s: %s\n", filename, strerror(errno) ); 00111 return MB_FILE_DOES_NOT_EXIST; 00112 } 00113 FileTokenizer tokens( file_ptr, readMeshIface ); 00114 00115 // Determine file format version 00116 const char* const start_tokens[] = { "$NOD", "$MeshFormat", 0 }; 00117 int format_version = tokens.match_token( start_tokens ); 00118 if (!format_version) 00119 return MB_FILE_DOES_NOT_EXIST; 00120 00121 // If version 2.0, read additional header info 00122 if (2 == format_version) 00123 { 00124 double version; 00125 if (!tokens.get_doubles( 1, &version )) 00126 return MB_FILE_WRITE_ERROR; 00127 00128 if (version != 2.0 && version != 2.1) { 00129 readMeshIface->report_error( "%s: unknown format version: %f\n", filename, version ); 00130 return MB_FILE_DOES_NOT_EXIST; 00131 } 00132 00133 int file_format; 00134 if (!tokens.get_integers( 1, &file_format ) || 00135 !tokens.get_integers( 1, &data_size ) || 00136 !tokens.match_token( "$EndMeshFormat" ) || 00137 !tokens.match_token( "$Nodes") ) 00138 return MB_FILE_WRITE_ERROR; 00139 } 00140 00141 // read number of nodes 00142 long num_nodes; 00143 if (!tokens.get_long_ints( 1, &num_nodes )) 00144 return MB_FILE_WRITE_ERROR; 00145 00146 // allocate nodes 00147 std::vector<double*> coord_arrays; 00148 EntityHandle handle = 0; 00149 result = readMeshIface->get_node_coords( 3, num_nodes, MB_START_ID, 00150 handle, coord_arrays ); 00151 if (MB_SUCCESS != result) 00152 return result; 00153 00154 // read nodes 00155 double *x = coord_arrays[0], 00156 *y = coord_arrays[1], 00157 *z = coord_arrays[2]; 00158 for( long i = 0; i < num_nodes; ++i, ++handle ) 00159 { 00160 long id; 00161 if (!tokens.get_long_ints( 1, &id ) || 00162 !tokens.get_doubles( 1, x++ ) || 00163 !tokens.get_doubles( 1, y++ ) || 00164 !tokens.get_doubles( 1, z++ )) 00165 return MB_FILE_WRITE_ERROR; 00166 00167 if (!node_id_map.insert( std::pair<long,EntityHandle>( id, handle ) ).second) 00168 { 00169 readMeshIface->report_error( "Dulicate node ID at line %d\n", 00170 tokens.line_number() ); 00171 return MB_FILE_WRITE_ERROR; 00172 } 00173 } 00174 00175 // create reverse map from handle to id 00176 std::vector<int> ids( num_nodes ); 00177 std::vector<int>::iterator id_iter = ids.begin(); 00178 std::vector<EntityHandle> handles( num_nodes ); 00179 std::vector<EntityHandle>::iterator h_iter = handles.begin(); 00180 for (std::map<long,EntityHandle>::iterator i = node_id_map.begin(); 00181 i != node_id_map.end(); ++i, ++id_iter, ++h_iter) 00182 { 00183 *id_iter = i->first; 00184 * h_iter = i->second; 00185 } 00186 // store IDs in tags 00187 result = mdbImpl->tag_set_data( globalId, &handles[0], num_nodes, &ids[0] ); 00188 if (MB_SUCCESS != result) 00189 return result; 00190 if (file_id_tag) { 00191 result = mdbImpl->tag_set_data( *file_id_tag, &handles[0], num_nodes, &ids[0] ); 00192 if (MB_SUCCESS != result) 00193 return result; 00194 } 00195 ids.clear(); 00196 handles.clear(); 00197 00198 // get tokens signifying end of node data and start of elements 00199 if (!tokens.match_token( format_version == 1 ? "$ENDNOD" : "$EndNodes" ) || 00200 !tokens.match_token( format_version == 1 ? "$ELM" : "$Elements" )) 00201 return MB_FILE_WRITE_ERROR; 00202 00203 // get element count 00204 long num_elem; 00205 if (!tokens.get_long_ints( 1, &num_elem )) 00206 return MB_FILE_WRITE_ERROR; 00207 00208 // lists of data accumulated for elements 00209 std::vector<EntityHandle> connectivity; 00210 std::vector<int> mat_set_list, geom_set_list, part_set_list, id_list; 00211 // temporary, per-element data 00212 std::vector<int> int_data(5), tag_data(2); 00213 std::vector<long> tmp_conn; 00214 int curr_elem_type = -1; 00215 for (long i = 0; i < num_elem; ++i) 00216 { 00217 // Read element description 00218 // File format 1.0 00219 if (1 == format_version) 00220 { 00221 if (!tokens.get_integers( 5, &int_data[0] )) 00222 return MB_FILE_WRITE_ERROR; 00223 tag_data[0] = int_data[2]; 00224 tag_data[1] = int_data[3]; 00225 if ((unsigned)tag_data[1] < GmshUtil::numGmshElemType && 00226 GmshUtil::gmshElemTypes[tag_data[1]].num_nodes != (unsigned)int_data[4]) { 00227 readMeshIface->report_error( "Invalid node count for element type at line %d\n", 00228 tokens.line_number() ); 00229 return MB_FILE_WRITE_ERROR; 00230 } 00231 } 00232 // File format 2.0 00233 else 00234 { 00235 if (!tokens.get_integers( 3, &int_data[0] )) 00236 return MB_FILE_WRITE_ERROR; 00237 tag_data.resize( int_data[2] ); 00238 if (!tokens.get_integers( tag_data.size(), &tag_data[0] )) 00239 return MB_FILE_WRITE_ERROR; 00240 } 00241 00242 // If a list of material sets was specified in the 00243 // argument list, skip any elements for which the 00244 // material set is not specified or is not in the 00245 // passed list. 00246 if (!blocks.empty() && (tag_data.empty() || 00247 blocks.find( tag_data[0] ) != blocks.end())) 00248 continue; 00249 00250 00251 // If the next element is not the same type as the last one, 00252 // create a sequence for the block of elements we've read 00253 // to this point (all of the same type), and clear accumulated 00254 // data. 00255 if (int_data[1] != curr_elem_type) 00256 { 00257 if (!id_list.empty()) // first iteration 00258 { 00259 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], 00260 id_list, 00261 mat_set_list, 00262 geom_set_list, 00263 part_set_list, 00264 connectivity, 00265 file_id_tag ) ; 00266 if (MB_SUCCESS != result) 00267 return result; 00268 } 00269 00270 id_list.clear(); 00271 mat_set_list.clear(); 00272 geom_set_list.clear(); 00273 part_set_list.clear(); 00274 connectivity.clear(); 00275 curr_elem_type = int_data[1]; 00276 if ((unsigned)curr_elem_type >= GmshUtil::numGmshElemType || 00277 GmshUtil::gmshElemTypes[curr_elem_type].mb_type == MBMAXTYPE) 00278 { 00279 readMeshIface->report_error( "Unsupported element type %d at line %d\n", 00280 curr_elem_type, tokens.line_number() ); 00281 return MB_FILE_WRITE_ERROR; 00282 } 00283 tmp_conn.resize( GmshUtil::gmshElemTypes[curr_elem_type].num_nodes ); 00284 } 00285 00286 // Store data from element description 00287 id_list.push_back( int_data[0] ); 00288 part_set_list.push_back( tag_data.size() > 2 ? tag_data[2] : 0 ); 00289 geom_set_list.push_back( tag_data.size() > 1 ? tag_data[1] : 0 ); 00290 mat_set_list.push_back( tag_data.size() > 0 ? tag_data[0] : 0 ); 00291 00292 // Get element connectivity 00293 if (!tokens.get_long_ints( tmp_conn.size(), &tmp_conn[0] )) 00294 return MB_FILE_WRITE_ERROR; 00295 00296 // Convert conectivity from IDs to handles 00297 for (unsigned j = 0; j < tmp_conn.size(); ++j) 00298 { 00299 std::map<long,EntityHandle>::iterator k = node_id_map.find( tmp_conn[j] ); 00300 if (k == node_id_map.end()) { 00301 readMeshIface->report_error( "Invalid node ID at line %d\n", 00302 tokens.line_number() ); 00303 return MB_FILE_WRITE_ERROR; 00304 } 00305 connectivity.push_back( k->second ); 00306 } 00307 } // for (num_nodes) 00308 00309 // Create entity sequence for last element(s). 00310 if (!id_list.empty()) 00311 { 00312 result = create_elements( GmshUtil::gmshElemTypes[curr_elem_type], 00313 id_list, 00314 mat_set_list, 00315 geom_set_list, 00316 part_set_list, 00317 connectivity, 00318 file_id_tag ) ; 00319 if (MB_SUCCESS != result) 00320 return result; 00321 } 00322 00323 // Construct parent-child relations for geometric sets. 00324 // Note: At the time this comment was written, the following 00325 // function was not impelemented. 00326 result = create_geometric_topology(); 00327 geomSets.clear(); 00328 return result; 00329 } 00330 00332 ErrorCode ReadGmsh::create_elements( const GmshElemType& type, 00333 const std::vector<int>& elem_ids, 00334 const std::vector<int>& matl_ids, 00335 const std::vector<int>& geom_ids, 00336 const std::vector<int>& prtn_ids, 00337 const std::vector<EntityHandle>& connectivity, 00338 const Tag* file_id_tag ) 00339 { 00340 ErrorCode result; 00341 00342 // Make sure input is consistent 00343 const unsigned long num_elem = elem_ids.size(); 00344 const int node_per_elem = type.num_nodes; 00345 if (matl_ids.size() != num_elem || 00346 geom_ids.size() != num_elem || 00347 prtn_ids.size() != num_elem || 00348 connectivity.size() != num_elem*node_per_elem) 00349 return MB_FAILURE; 00350 00351 // Create the element sequence 00352 EntityHandle handle = 0; 00353 EntityHandle* conn_array; 00354 result = readMeshIface->get_element_connect( num_elem, node_per_elem, type.mb_type, 00355 MB_START_ID, 00356 handle, conn_array ); 00357 if (MB_SUCCESS != result) 00358 return result; 00359 00360 // Copy passed element connectivity into entity sequence data. 00361 if (type.node_order) 00362 { 00363 for (unsigned long i = 0; i < num_elem; ++i) 00364 for (int j = 0; j < node_per_elem; ++j) 00365 conn_array[i*node_per_elem+type.node_order[j]] = connectivity[i*node_per_elem+j]; 00366 } 00367 else 00368 { 00369 memcpy( conn_array, &connectivity[0], connectivity.size() * sizeof(EntityHandle) ); 00370 } 00371 00372 // notify MOAB of the new elements 00373 result = readMeshIface->update_adjacencies(handle, num_elem, node_per_elem, conn_array); 00374 if (MB_SUCCESS != result) return result; 00375 00376 // Store element IDs 00377 Range elements( handle, handle + num_elem - 1 ); 00378 result = mdbImpl->tag_set_data( globalId, elements, &elem_ids[0] ); 00379 if (MB_SUCCESS != result) 00380 return result; 00381 if (file_id_tag) { 00382 result = mdbImpl->tag_set_data( *file_id_tag, elements, &elem_ids[0] ); 00383 if (MB_SUCCESS != result) 00384 return result; 00385 } 00386 00387 // Add elements to material sets 00388 result = create_sets( type.mb_type, elements, matl_ids, 0 ); 00389 if (MB_SUCCESS != result) 00390 return result; 00391 // Add elements to geometric sets 00392 result = create_sets( type.mb_type, elements, geom_ids, 1 ); 00393 if (MB_SUCCESS != result) 00394 return result; 00395 // Add elements to parallel partitions 00396 result = create_sets( type.mb_type, elements, prtn_ids, 2 ); 00397 if (MB_SUCCESS != result) 00398 return result; 00399 00400 return MB_SUCCESS; 00401 } 00402 00404 ErrorCode ReadGmsh::create_sets( EntityType type, 00405 const Range& elements, 00406 const std::vector<int>& set_ids, 00407 int set_type ) 00408 { 00409 ErrorCode result; 00410 00411 // Get a unque list of set IDs 00412 std::set<int> ids; 00413 for (std::vector<int>::const_iterator i = set_ids.begin(); i != set_ids.end(); ++i) 00414 ids.insert( *i ); 00415 00416 // No Sets? 00417 if (ids.empty() || (ids.size() == 1 && *ids.begin() == 0)) 00418 return MB_SUCCESS; // no sets (all ids are zero) 00419 00420 00421 // Get/create tag handles 00422 int num_tags; 00423 Tag tag_handles[2]; 00424 int tag_val; 00425 const void* tag_values[2] = { &tag_val, NULL }; 00426 00427 switch( set_type ) 00428 { 00429 default: 00430 return MB_FAILURE; 00431 case 0: 00432 case 2: 00433 { 00434 const char* name = set_type ? PARALLEL_PARTITION_TAG_NAME : MATERIAL_SET_TAG_NAME; 00435 result = mdbImpl->tag_get_handle( name, 1, MB_TYPE_INTEGER, tag_handles[0], 00436 MB_TAG_SPARSE|MB_TAG_CREAT ); 00437 if (MB_SUCCESS != result) 00438 return result; 00439 num_tags = 1; 00440 break; 00441 } 00442 case 1: 00443 { 00444 result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, 00445 tag_handles[1], MB_TAG_SPARSE|MB_TAG_CREAT ); 00446 if (MB_SUCCESS != result) 00447 return result; 00448 tag_values[1] = NULL; 00449 tag_handles[0]= globalId; 00450 num_tags = 2; 00451 break; 00452 } 00453 } // switch 00454 00455 // For each unique set ID... 00456 for (std::set<int>::iterator i = ids.begin(); i != ids.end(); ++i) 00457 { 00458 // Skip "null" set ID 00459 if (*i == 0) 00460 continue; 00461 00462 // Get all entities with the current set ID 00463 Range entities, sets; 00464 std::vector<int>::const_iterator j = set_ids.begin(); 00465 for (Range::iterator k = elements.begin(); k != elements.end(); ++j, ++k) 00466 if (*i == *j) 00467 entities.insert( *k ); 00468 00469 // Get set by ID 00470 tag_val = *i; 00471 result = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, 00472 tag_handles, tag_values, num_tags, 00473 sets ); 00474 if (MB_SUCCESS != result && MB_ENTITY_NOT_FOUND != result) 00475 return result; 00476 00477 // Don't use existing geometry sets (from some other file) 00478 if (1 == set_type) // geometry 00479 sets = intersect( sets, geomSets ); 00480 00481 // Get set handle 00482 EntityHandle set; 00483 // If no sets with ID, create one 00484 if (sets.empty()) 00485 { 00486 result = mdbImpl->create_meshset( MESHSET_SET, set ); 00487 if (MB_SUCCESS != result) 00488 return result; 00489 00490 result = mdbImpl->tag_set_data( tag_handles[0], &set, 1, &*i ); 00491 if (MB_SUCCESS != result) 00492 return result; 00493 00494 if (1 == set_type) // geometry 00495 { 00496 int dim = CN::Dimension( type ); 00497 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim ); 00498 if (MB_SUCCESS != result) 00499 return result; 00500 geomSets.insert( set ); 00501 } 00502 } 00503 else 00504 { 00505 set = *sets.begin(); 00506 if (1 == set_type) // geometry 00507 { 00508 int dim = CN::Dimension( type ); 00509 // Get dimension of set 00510 int dim2; 00511 result = mdbImpl->tag_get_data( tag_handles[1], &set, 1, &dim2 ); 00512 if (MB_SUCCESS != result) 00513 return result; 00514 // If we're putting geometry of a higher dimension into the 00515 // set, increase the dimension of the set. 00516 if (dim > dim2) { 00517 result = mdbImpl->tag_set_data( tag_handles[1], &set, 1, &dim ); 00518 if (MB_SUCCESS != result) 00519 return result; 00520 } 00521 } 00522 } 00523 00524 // Put the mesh entities into the set 00525 result = mdbImpl->add_entities( set, entities ); 00526 if (MB_SUCCESS != result) 00527 return result; 00528 } // for (ids) 00529 00530 return MB_SUCCESS; 00531 } 00532 00533 00537 ErrorCode ReadGmsh::create_geometric_topology() 00538 { 00539 if (geomSets.empty()) 00540 return MB_SUCCESS; 00541 00542 // not implemented yet 00543 geomSets.clear(); 00544 return MB_SUCCESS; 00545 } 00546 00547 } // namespace moab