moab
|
00001 00022 #include "ReadSms.hpp" 00023 #include "FileTokenizer.hpp" // for file tokenizer 00024 #include "Internals.hpp" 00025 #include "moab/Interface.hpp" 00026 #include "moab/ReadUtilIface.hpp" 00027 #include "moab/Range.hpp" 00028 #include "MBTagConventions.hpp" 00029 #include "MBParallelConventions.h" 00030 #include "moab/CN.hpp" 00031 00032 #include <errno.h> 00033 #include <string.h> 00034 #include <map> 00035 #include <set> 00036 #include <iostream> 00037 00038 #define CHECK(a) if (MB_SUCCESS != result) { \ 00039 std::cerr << a << std::endl; \ 00040 return result; \ 00041 } 00042 00043 #define CHECKN(a) if (n != (a)) return MB_FILE_WRITE_ERROR 00044 00045 namespace moab { 00046 00047 ReaderIface* ReadSms::factory( Interface* iface ) 00048 { return new ReadSms(iface); } 00049 00050 ReadSms::ReadSms(Interface* impl) 00051 : mdbImpl(impl) 00052 { 00053 mdbImpl->query_interface(readMeshIface); 00054 } 00055 00056 ReadSms::~ReadSms() 00057 { 00058 if (readMeshIface) { 00059 mdbImpl->release_interface(readMeshIface); 00060 readMeshIface = 0; 00061 } 00062 } 00063 00064 ErrorCode ReadSms::read_tag_values(const char* /* file_name */, 00065 const char* /* tag_name */, 00066 const FileOptions& /* opts */, 00067 std::vector<int>& /* tag_values_out */, 00068 const SubsetList* /* subset_list */ ) 00069 { 00070 return MB_NOT_IMPLEMENTED; 00071 } 00072 00073 00074 ErrorCode ReadSms::load_file( const char* filename, 00075 const EntityHandle* , 00076 const FileOptions& , 00077 const ReaderIface::SubsetList* subset_list, 00078 const Tag* file_id_tag ) 00079 { 00080 if (subset_list) { 00081 readMeshIface->report_error( "Reading subset of files not supported for Sms." ); 00082 return MB_UNSUPPORTED_OPERATION; 00083 } 00084 00085 setId = 1; 00086 00087 // Open file 00088 FILE* file_ptr = fopen( filename, "r" ); 00089 if (!file_ptr) 00090 { 00091 readMeshIface->report_error( "%s: %s\n", filename, strerror(errno) ); 00092 return MB_FILE_DOES_NOT_EXIST; 00093 } 00094 00095 const ErrorCode result = load_file_impl( file_ptr, file_id_tag ); 00096 fclose( file_ptr ); 00097 00098 return result; 00099 } 00100 00101 ErrorCode ReadSms::load_file_impl( FILE* file_ptr, const Tag* file_id_tag ) 00102 { 00103 bool warned = false; 00104 double dum_params[] = {0.0, 0.0, 0.0}; 00105 00106 int zero = 0; 00107 ErrorCode result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, 00108 globalId, MB_TAG_DENSE|MB_TAG_CREAT, &zero); 00109 CHECK("Failed to create gid tag."); 00110 00111 result = mdbImpl->tag_get_handle("PARAMETER_COORDS", 3, MB_TYPE_DOUBLE, 00112 paramCoords, MB_TAG_DENSE|MB_TAG_CREAT ); 00113 CHECK("Failed to create param coords tag."); 00114 00115 int negone = -1; 00116 result = mdbImpl->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, 00117 geomDimension, MB_TAG_SPARSE|MB_TAG_CREAT, &negone); 00118 CHECK("Failed to create geom dim tag."); 00119 00120 int n; 00121 char line[256]; 00122 int file_type; 00123 n = fscanf(file_ptr, "%s %d", line, &file_type); 00124 CHECKN(2); 00125 00126 if (3 == file_type) { 00127 result = read_parallel_info(file_ptr); 00128 CHECK("Failed to read parallel info."); 00129 } 00130 00131 int nregions, nfaces, nedges, nvertices, npoints; 00132 n = fscanf(file_ptr, "%d %d %d %d %d", &nregions, &nfaces, &nedges, 00133 &nvertices, &npoints); 00134 CHECKN(5); 00135 if (nregions < 0 || nfaces < 0 || nedges < 0 || nvertices < 0 || npoints < 0) 00136 return MB_FILE_WRITE_ERROR; 00137 00138 // create the vertices 00139 std::vector<double*> coord_arrays; 00140 EntityHandle vstart = 0; 00141 result = readMeshIface->get_node_coords( 3, nvertices, MB_START_ID, 00142 vstart, coord_arrays ); 00143 CHECK("Failed to get node arrays."); 00144 if (MB_SUCCESS != result) 00145 return result; 00146 00147 if (file_id_tag) { 00148 result = add_entities( vstart, nvertices, file_id_tag ); 00149 if (MB_SUCCESS != result) 00150 return result; 00151 } 00152 00153 EntityHandle this_gent, new_handle; 00154 std::vector<EntityHandle> gentities[4]; 00155 int gent_id, dum_int; 00156 int gent_type, num_connections; 00157 00158 for(int i = 0; i < nvertices; i++) 00159 { 00160 n = fscanf(file_ptr, "%d", &gent_id); 00161 CHECKN(1); 00162 if (!gent_id) continue; 00163 00164 n = fscanf(file_ptr,"%d %d %lf %lf %lf", &gent_type, &num_connections, 00165 coord_arrays[0]+i, coord_arrays[1]+i, coord_arrays[2]+i); 00166 CHECKN(5); 00167 00168 result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag ); 00169 if (MB_SUCCESS != result) 00170 return result; 00171 00172 new_handle = vstart + i; 00173 result = mdbImpl->add_entities(this_gent, &new_handle, 1); 00174 CHECK("Adding vertex to geom set failed."); 00175 if (MB_SUCCESS != result) return result; 00176 00177 switch(gent_type) 00178 { 00179 case 1: 00180 n = fscanf(file_ptr, "%le", dum_params); 00181 CHECKN(1); 00182 result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params); 00183 CHECK("Failed to set param coords tag for vertex."); 00184 if (MB_SUCCESS != result) return result; 00185 break; 00186 case 2: 00187 n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int); 00188 CHECKN(3); 00189 dum_params[2] = dum_int; 00190 result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params); 00191 CHECK("Failed to set param coords tag for vertex."); 00192 if (MB_SUCCESS != result) return result; 00193 break; 00194 00195 default: break; 00196 } 00197 } // end of reading vertices 00198 00199 // ******************************* 00200 // Read Edges 00201 // ******************************* 00202 00203 int vert1, vert2, num_pts; 00204 std::vector<EntityHandle> everts(2); 00205 EntityHandle estart, *connect; 00206 result = readMeshIface->get_element_connect(nedges, 2, MBEDGE, 1, estart, connect); 00207 CHECK("Failed to create array of edges."); 00208 if (MB_SUCCESS != result) return result; 00209 00210 if (file_id_tag) { 00211 result = add_entities( estart, nedges, file_id_tag ); 00212 if (MB_SUCCESS != result) 00213 return result; 00214 } 00215 00216 for(int i = 0; i < nedges; i++) 00217 { 00218 n = fscanf(file_ptr,"%d",&gent_id); 00219 CHECKN(1); 00220 if (!gent_id) continue; 00221 00222 n = fscanf(file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, 00223 &num_connections, &num_pts); 00224 CHECKN(5); 00225 if (vert1 < 1 || vert1 > nvertices) 00226 return MB_FILE_WRITE_ERROR; 00227 if (vert2 < 1 || vert2 > nvertices) 00228 return MB_FILE_WRITE_ERROR; 00229 00230 connect[0] = vstart + vert1 - 1; 00231 connect[1] = vstart + vert2 - 1; 00232 if (num_pts > 1 && !warned) { 00233 std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl; 00234 warned = true; 00235 } 00236 00237 result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag); 00238 CHECK("Problem getting geom set for edge."); 00239 if (MB_SUCCESS != result) 00240 return result; 00241 00242 new_handle = estart + i; 00243 result = mdbImpl->add_entities(this_gent, &new_handle, 1); 00244 CHECK("Failed to add edge to geom set."); 00245 if (MB_SUCCESS != result) return result; 00246 00247 connect += 2; 00248 00249 for(int j = 0; j < num_pts; j++) { 00250 switch(gent_type) { 00251 case 1: 00252 n = fscanf(file_ptr, "%le", dum_params); 00253 CHECKN(1); 00254 result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params); 00255 CHECK("Failed to set param coords tag for edge."); 00256 if (MB_SUCCESS != result) return result; 00257 break; 00258 case 2: 00259 n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int); 00260 CHECKN(3); 00261 dum_params[2] = dum_int; 00262 result = mdbImpl->tag_set_data(paramCoords, &new_handle, 1, dum_params); 00263 CHECK("Failed to set param coords tag for edge."); 00264 if (MB_SUCCESS != result) return result; 00265 break; 00266 default: 00267 break; 00268 } 00269 } 00270 00271 } // end of reading edges 00272 00273 // ******************************* 00274 // Read Faces 00275 // ******************************* 00276 std::vector<EntityHandle> bound_ents, bound_verts, new_faces; 00277 int bound_id; 00278 Range shverts; 00279 new_faces.resize(nfaces); 00280 int num_bounding; 00281 00282 for(int i = 0; i < nfaces; i++) 00283 { 00284 n = fscanf(file_ptr, "%d", &gent_id); 00285 CHECKN(1); 00286 if(!gent_id) continue; 00287 00288 n = fscanf(file_ptr,"%d %d", &gent_type, &num_bounding); 00289 CHECKN(2); 00290 00291 result = get_set(gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag); 00292 CHECK("Problem getting geom set for face."); 00293 if (MB_SUCCESS != result) 00294 return result; 00295 00296 bound_ents.resize(num_bounding+1); 00297 bound_verts.resize(num_bounding); 00298 for(int j = 0; j < num_bounding; j++) { 00299 n = fscanf(file_ptr, "%d ", &bound_id); 00300 CHECKN(1); 00301 if (0 > bound_id) bound_id = abs(bound_id); 00302 assert(0 < bound_id && bound_id <= nedges); 00303 if (bound_id < 1 || bound_id > nedges) 00304 return MB_FILE_WRITE_ERROR; 00305 bound_ents[j] = estart + abs(bound_id) - 1; 00306 } 00307 00308 // convert edge-based model to vertex-based one 00309 for (int j = 0; j < num_bounding; j++) { 00310 if (j == num_bounding-1) bound_ents[j+1] = bound_ents[0]; 00311 result = mdbImpl->get_adjacencies(&bound_ents[j], 2, 0, false, shverts); 00312 CHECK("Failed to get vertices bounding edge."); 00313 if (MB_SUCCESS != result) return result; 00314 assert(shverts.size() == 1); 00315 bound_verts[j] = *shverts.begin(); 00316 shverts.clear(); 00317 } 00318 00319 result = mdbImpl->create_element((EntityType)(MBTRI+num_bounding-3), 00320 &bound_verts[0], bound_verts.size(), 00321 new_faces[i]); 00322 CHECK("Failed to create edge."); 00323 if (MB_SUCCESS != result) return result; 00324 00325 result = mdbImpl->add_entities(this_gent, &new_faces[i], 1); 00326 CHECK("Failed to add edge to geom set."); 00327 if (MB_SUCCESS != result) return result; 00328 00329 int num_read = fscanf(file_ptr, "%d", &num_pts); 00330 if(!num_pts || !num_read) continue; 00331 00332 for(int j = 0; j < num_pts; j++) { 00333 switch(gent_type) { 00334 case 1: 00335 n = fscanf(file_ptr, "%le", dum_params); CHECKN(1); 00336 result = mdbImpl->tag_set_data(paramCoords, &new_faces[i], 1, dum_params); 00337 CHECK("Failed to set param coords tag for face."); 00338 if (MB_SUCCESS != result) return result; 00339 break; 00340 case 2: 00341 n = fscanf(file_ptr, "%le %le %d", dum_params, dum_params+1, &dum_int); 00342 CHECKN(3); 00343 dum_params[2] = dum_int; 00344 result = mdbImpl->tag_set_data(paramCoords, &new_faces[i], 1, dum_params); 00345 CHECK("Failed to set param coords tag for face."); 00346 if (MB_SUCCESS != result) return result; 00347 break; 00348 default: 00349 break; 00350 } 00351 } 00352 00353 } // end of reading faces 00354 00355 if (file_id_tag) { 00356 result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 ); 00357 if (MB_SUCCESS != result) 00358 return result; 00359 } 00360 00361 00362 // ******************************* 00363 // Read Regions 00364 // ******************************* 00365 int sense[MAX_SUB_ENTITIES]; 00366 bound_verts.resize(MAX_SUB_ENTITIES); 00367 00368 std::vector<EntityHandle> regions; 00369 if (file_id_tag) 00370 regions.resize( nregions ); 00371 for(int i = 0; i < nregions; i++) 00372 { 00373 n = fscanf(file_ptr, "%d", &gent_id); CHECKN(1); 00374 if (!gent_id) continue; 00375 result = get_set(gentities, 3, gent_id, geomDimension, this_gent, file_id_tag); 00376 CHECK("Couldn't get geom set for region."); 00377 if (MB_SUCCESS != result) 00378 return result; 00379 n = fscanf(file_ptr, "%d", &num_bounding); CHECKN(1); 00380 bound_ents.resize(num_bounding); 00381 for(int j = 0; j < num_bounding; j++) { 00382 n = fscanf(file_ptr, "%d ", &bound_id); CHECKN(1); 00383 assert(abs(bound_id) < (int)new_faces.size()+1 && bound_id); 00384 if (!bound_id || abs(bound_id) > nfaces) 00385 return MB_FILE_WRITE_ERROR; 00386 sense[j] = (bound_id < 0) ? -1 : 1; 00387 bound_ents[j] = new_faces[abs(bound_id)-1]; 00388 } 00389 00390 EntityType etype; 00391 result = readMeshIface->get_ordered_vertices(&bound_ents[0], sense, 00392 num_bounding, 00393 3, &bound_verts[0], etype); 00394 CHECK("Failed in get_ordered_vertices."); 00395 00396 // make the element 00397 result = mdbImpl->create_element(etype, &bound_verts[0], 00398 CN::VerticesPerEntity(etype), new_handle); 00399 CHECK("Failed to create region."); 00400 if (MB_SUCCESS != result) return result; 00401 00402 result = mdbImpl->add_entities(this_gent, &new_handle, 1); 00403 CHECK("Failed to add region to geom set."); 00404 if (MB_SUCCESS != result) return result; 00405 00406 if (file_id_tag) 00407 regions[i] = new_handle; 00408 00409 n = fscanf(file_ptr, "%d ", &dum_int); CHECKN(1); 00410 00411 } // end of reading regions 00412 00413 if (file_id_tag) { 00414 result = readMeshIface->assign_ids( *file_id_tag, ®ions[0], regions.size(), 1 ); 00415 if (MB_SUCCESS != result) 00416 return result; 00417 } 00418 00419 return MB_SUCCESS; 00420 } 00421 00422 ErrorCode ReadSms::get_set(std::vector<EntityHandle> *sets, 00423 int set_dim, int set_id, 00424 Tag dim_tag, 00425 EntityHandle &this_set, 00426 const Tag* file_id_tag) 00427 { 00428 ErrorCode result = MB_SUCCESS; 00429 00430 if (set_dim < 0 || set_dim > 3) 00431 return MB_FILE_WRITE_ERROR; 00432 00433 if ((int)sets[set_dim].size() <= set_id || 00434 !sets[set_dim][set_id]) { 00435 if ((int)sets[set_dim].size() <= set_id) 00436 sets[set_dim].resize(set_id+1, 0); 00437 00438 if (!sets[set_dim][set_id]) { 00439 result = mdbImpl->create_meshset(MESHSET_SET, 00440 sets[set_dim][set_id]); 00441 if (MB_SUCCESS != result) return result; 00442 result = mdbImpl->tag_set_data(globalId, 00443 &sets[set_dim][set_id], 1, 00444 &set_id); 00445 if (MB_SUCCESS != result) return result; 00446 result = mdbImpl->tag_set_data(dim_tag, 00447 &sets[set_dim][set_id], 1, 00448 &set_dim); 00449 if (MB_SUCCESS != result) return result; 00450 00451 if (file_id_tag) { 00452 result = mdbImpl->tag_set_data(*file_id_tag, 00453 &sets[set_dim][set_id], 1, 00454 &setId); 00455 ++setId; 00456 } 00457 } 00458 } 00459 00460 this_set = sets[set_dim][set_id]; 00461 00462 return result; 00463 } 00464 00465 ErrorCode ReadSms::read_parallel_info(FILE *file_ptr) 00466 { 00467 // ErrorCode result; 00468 00469 // read partition info 00470 int nparts, part_id, num_ifaces, num_corner_ents; 00471 int num_read = fscanf(file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents); 00472 if (!num_read) return MB_FAILURE; 00473 00474 // read interfaces 00475 int iface_id, iface_dim, iface_own, num_iface_corners; 00476 // EntityHandle this_iface; 00477 std::vector<int> *iface_corners; 00478 for (int i = 0; i < num_ifaces; i++) { 00479 num_read = fscanf(file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own, 00480 &num_iface_corners); 00481 if (!num_read) return MB_FAILURE; 00482 00483 // result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface); 00484 // CHECK("Failed to make iface set."); 00485 00486 // read the corner ids and store them on the set for now 00487 iface_corners = new std::vector<int>(num_iface_corners); 00488 for (int j = 0; j < num_iface_corners; j++) { 00489 num_read = fscanf(file_ptr, "%d", &(*iface_corners)[j]); 00490 if (!num_read) return MB_FAILURE; 00491 } 00492 00493 // result = tag_set_data(ifaceCornerTag, &this_iface, 1, 00494 // &iface_corners); 00495 // CHECK("Failed to set iface corner tag."); 00496 } 00497 00498 // interface data has been read 00499 return MB_SUCCESS; 00500 } 00501 00502 ErrorCode ReadSms::add_entities( EntityHandle start, 00503 EntityHandle count, 00504 const Tag* file_id_tag ) 00505 { 00506 if (!count || !file_id_tag) 00507 return MB_FAILURE; 00508 00509 Range range; 00510 range.insert( start, start + count - 1 ); 00511 return readMeshIface->assign_ids( *file_id_tag, range, 1 ); 00512 } 00513 00514 } // namespace moab