MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in this software. 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 2.1 of the License, or (at your option) any later version. 00013 * 00014 */ 00015 00016 /** 00017 * \class ReadSms 00018 * \brief Sms (http://www.geuz.org/sms) file reader 00019 * \author Jason Kraftcheck 00020 */ 00021 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 ) \ 00039 if( MB_SUCCESS != result ) \ 00040 { \ 00041 std::cerr << a << std::endl; \ 00042 return result; \ 00043 } 00044 00045 #define CHECKN( a ) \ 00046 if( n != ( a ) ) return MB_FILE_WRITE_ERROR 00047 00048 namespace moab 00049 { 00050 00051 ReaderIface* ReadSms::factory( Interface* iface ) 00052 { 00053 return new ReadSms( iface ); 00054 } 00055 00056 ReadSms::ReadSms( Interface* impl ) : mdbImpl( impl ), globalId( 0 ), paramCoords( 0 ), geomDimension( 0 ), setId( 0 ) 00057 { 00058 mdbImpl->query_interface( readMeshIface ); 00059 } 00060 00061 ReadSms::~ReadSms() 00062 { 00063 if( readMeshIface ) 00064 { 00065 mdbImpl->release_interface( readMeshIface ); 00066 readMeshIface = 0; 00067 } 00068 } 00069 00070 ErrorCode ReadSms::read_tag_values( const char* /* file_name */, const char* /* tag_name */, 00071 const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */, 00072 const SubsetList* /* subset_list */ ) 00073 { 00074 return MB_NOT_IMPLEMENTED; 00075 } 00076 00077 ErrorCode ReadSms::load_file( const char* filename, const EntityHandle* /* file_set */, const FileOptions& /* opts */, 00078 const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag ) 00079 { 00080 if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for Sms" ); } 00081 00082 setId = 1; 00083 00084 // Open file 00085 FILE* file_ptr = fopen( filename, "r" ); 00086 if( !file_ptr ) { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) ); } 00087 00088 const ErrorCode result = load_file_impl( file_ptr, file_id_tag ); 00089 fclose( file_ptr ); 00090 00091 return result; 00092 } 00093 00094 ErrorCode ReadSms::load_file_impl( FILE* file_ptr, const Tag* file_id_tag ) 00095 { 00096 bool warned = false; 00097 double dum_params[] = { 0.0, 0.0, 0.0 }; 00098 00099 globalId = mdbImpl->globalId_tag(); 00100 00101 ErrorCode result = 00102 mdbImpl->tag_get_handle( "PARAMETER_COORDS", 3, MB_TYPE_DOUBLE, paramCoords, MB_TAG_DENSE | MB_TAG_CREAT ); 00103 CHECK( "Failed to create param coords tag." ); 00104 00105 int negone = -1; 00106 result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomDimension, 00107 MB_TAG_SPARSE | MB_TAG_CREAT, &negone ); 00108 CHECK( "Failed to create geom dim tag." ); 00109 00110 int n; 00111 char line[256], all_line[1024]; 00112 int file_type; 00113 00114 if( fgets( all_line, sizeof( all_line ), file_ptr ) == NULL ) { return MB_FAILURE; } 00115 if( sscanf( all_line, "%s %d", line, &file_type ) != 2 ) { return MB_FAILURE; } 00116 00117 if( 3 == file_type ) 00118 { 00119 result = read_parallel_info( file_ptr ); 00120 CHECK( "Failed to read parallel info." ); 00121 } 00122 00123 int nregions, nfaces, nedges, nvertices, npoints; 00124 n = fscanf( file_ptr, "%d %d %d %d %d", &nregions, &nfaces, &nedges, &nvertices, &npoints ); 00125 CHECKN( 5 ); 00126 if( nregions < 0 || nfaces < 0 || nedges < 0 || nvertices < 0 || npoints < 0 ) return MB_FILE_WRITE_ERROR; 00127 00128 // Create the vertices 00129 std::vector< double* > coord_arrays; 00130 EntityHandle vstart = 0; 00131 result = readMeshIface->get_node_coords( 3, nvertices, MB_START_ID, vstart, coord_arrays ); 00132 CHECK( "Failed to get node arrays." ); 00133 00134 if( file_id_tag ) 00135 { 00136 result = add_entities( vstart, nvertices, file_id_tag );MB_CHK_ERR( result ); 00137 } 00138 00139 EntityHandle this_gent, new_handle; 00140 std::vector< EntityHandle > gentities[4]; 00141 int gent_id, dum_int; 00142 int gent_type, num_connections; 00143 00144 for( int i = 0; i < nvertices; i++ ) 00145 { 00146 n = fscanf( file_ptr, "%d", &gent_id ); 00147 CHECKN( 1 ); 00148 if( !gent_id ) continue; 00149 00150 n = fscanf( file_ptr, "%d %d %lf %lf %lf", &gent_type, &num_connections, coord_arrays[0] + i, 00151 coord_arrays[1] + i, coord_arrays[2] + i ); 00152 CHECKN( 5 ); 00153 00154 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );MB_CHK_ERR( result ); 00155 00156 new_handle = vstart + i; 00157 result = mdbImpl->add_entities( this_gent, &new_handle, 1 ); 00158 CHECK( "Adding vertex to geom set failed." ); 00159 00160 switch( gent_type ) 00161 { 00162 case 1: 00163 n = fscanf( file_ptr, "%le", dum_params ); 00164 CHECKN( 1 ); 00165 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params ); 00166 CHECK( "Failed to set param coords tag for vertex." ); 00167 break; 00168 case 2: 00169 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int ); 00170 CHECKN( 3 ); 00171 dum_params[2] = dum_int; 00172 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params ); 00173 CHECK( "Failed to set param coords tag for vertex." ); 00174 break; 00175 default: 00176 break; 00177 } 00178 } // End of reading vertices 00179 00180 // ******************************* 00181 // Read Edges 00182 // ******************************* 00183 00184 int vert1, vert2, num_pts; 00185 std::vector< EntityHandle > everts( 2 ); 00186 EntityHandle estart, *connect; 00187 result = readMeshIface->get_element_connect( nedges, 2, MBEDGE, 1, estart, connect ); 00188 CHECK( "Failed to create array of edges." ); 00189 00190 if( file_id_tag ) 00191 { 00192 result = add_entities( estart, nedges, file_id_tag ); 00193 if( MB_SUCCESS != result ) return result; 00194 } 00195 00196 for( int i = 0; i < nedges; i++ ) 00197 { 00198 n = fscanf( file_ptr, "%d", &gent_id ); 00199 CHECKN( 1 ); 00200 if( !gent_id ) continue; 00201 00202 n = fscanf( file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, &num_connections, &num_pts ); 00203 CHECKN( 5 ); 00204 if( vert1 < 1 || vert1 > nvertices ) return MB_FILE_WRITE_ERROR; 00205 if( vert2 < 1 || vert2 > nvertices ) return MB_FILE_WRITE_ERROR; 00206 00207 connect[0] = vstart + vert1 - 1; 00208 connect[1] = vstart + vert2 - 1; 00209 if( num_pts > 1 && !warned ) 00210 { 00211 std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl; 00212 warned = true; 00213 } 00214 00215 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag ); 00216 CHECK( "Problem getting geom set for edge." ); 00217 00218 new_handle = estart + i; 00219 result = mdbImpl->add_entities( this_gent, &new_handle, 1 ); 00220 CHECK( "Failed to add edge to geom set." ); 00221 00222 connect += 2; 00223 00224 for( int j = 0; j < num_pts; j++ ) 00225 { 00226 switch( gent_type ) 00227 { 00228 case 1: 00229 n = fscanf( file_ptr, "%le", dum_params ); 00230 CHECKN( 1 ); 00231 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params ); 00232 CHECK( "Failed to set param coords tag for edge." ); 00233 break; 00234 case 2: 00235 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int ); 00236 CHECKN( 3 ); 00237 dum_params[2] = dum_int; 00238 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params ); 00239 CHECK( "Failed to set param coords tag for edge." ); 00240 break; 00241 default: 00242 break; 00243 } 00244 } 00245 } // End of reading edges 00246 00247 // ******************************* 00248 // Read Faces 00249 // ******************************* 00250 std::vector< EntityHandle > bound_ents, bound_verts, new_faces; 00251 int bound_id; 00252 Range shverts; 00253 new_faces.resize( nfaces ); 00254 int num_bounding; 00255 00256 for( int i = 0; i < nfaces; i++ ) 00257 { 00258 n = fscanf( file_ptr, "%d", &gent_id ); 00259 CHECKN( 1 ); 00260 if( !gent_id ) continue; 00261 00262 n = fscanf( file_ptr, "%d %d", &gent_type, &num_bounding ); 00263 CHECKN( 2 ); 00264 00265 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag ); 00266 CHECK( "Problem getting geom set for face." ); 00267 00268 bound_ents.resize( num_bounding + 1 ); 00269 bound_verts.resize( num_bounding ); 00270 for( int j = 0; j < num_bounding; j++ ) 00271 { 00272 n = fscanf( file_ptr, "%d ", &bound_id ); 00273 CHECKN( 1 ); 00274 if( 0 > bound_id ) bound_id = abs( bound_id ); 00275 assert( 0 < bound_id && bound_id <= nedges ); 00276 if( bound_id < 1 || bound_id > nedges ) return MB_FILE_WRITE_ERROR; 00277 bound_ents[j] = estart + abs( bound_id ) - 1; 00278 } 00279 00280 // Convert edge-based model to vertex-based one 00281 for( int j = 0; j < num_bounding; j++ ) 00282 { 00283 if( j == num_bounding - 1 ) bound_ents[j + 1] = bound_ents[0]; 00284 result = mdbImpl->get_adjacencies( &bound_ents[j], 2, 0, false, shverts ); 00285 CHECK( "Failed to get vertices bounding edge." ); 00286 assert( shverts.size() == 1 ); 00287 bound_verts[j] = *shverts.begin(); 00288 shverts.clear(); 00289 } 00290 00291 result = mdbImpl->create_element( ( EntityType )( MBTRI + num_bounding - 3 ), &bound_verts[0], 00292 bound_verts.size(), new_faces[i] ); 00293 CHECK( "Failed to create edge." ); 00294 00295 result = mdbImpl->add_entities( this_gent, &new_faces[i], 1 ); 00296 CHECK( "Failed to add edge to geom set." ); 00297 00298 int num_read = fscanf( file_ptr, "%d", &num_pts ); 00299 if( !num_pts || !num_read ) continue; 00300 00301 for( int j = 0; j < num_pts; j++ ) 00302 { 00303 switch( gent_type ) 00304 { 00305 case 1: 00306 n = fscanf( file_ptr, "%le", dum_params ); 00307 CHECKN( 1 ); 00308 result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params ); 00309 CHECK( "Failed to set param coords tag for face." ); 00310 break; 00311 case 2: 00312 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int ); 00313 CHECKN( 3 ); 00314 dum_params[2] = dum_int; 00315 result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params ); 00316 CHECK( "Failed to set param coords tag for face." ); 00317 break; 00318 default: 00319 break; 00320 } 00321 } 00322 } // End of reading faces 00323 00324 if( file_id_tag ) 00325 { 00326 result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 ); 00327 if( MB_SUCCESS != result ) return result; 00328 } 00329 00330 // ******************************* 00331 // Read Regions 00332 // ******************************* 00333 int sense[MAX_SUB_ENTITIES]; 00334 bound_verts.resize( MAX_SUB_ENTITIES ); 00335 00336 std::vector< EntityHandle > regions; 00337 if( file_id_tag ) regions.resize( nregions ); 00338 for( int i = 0; i < nregions; i++ ) 00339 { 00340 n = fscanf( file_ptr, "%d", &gent_id ); 00341 CHECKN( 1 ); 00342 if( !gent_id ) continue; 00343 result = get_set( gentities, 3, gent_id, geomDimension, this_gent, file_id_tag ); 00344 CHECK( "Couldn't get geom set for region." ); 00345 n = fscanf( file_ptr, "%d", &num_bounding ); 00346 CHECKN( 1 ); 00347 bound_ents.resize( num_bounding ); 00348 for( int j = 0; j < num_bounding; j++ ) 00349 { 00350 n = fscanf( file_ptr, "%d ", &bound_id ); 00351 CHECKN( 1 ); 00352 assert( abs( bound_id ) < (int)new_faces.size() + 1 && bound_id ); 00353 if( !bound_id || abs( bound_id ) > nfaces ) return MB_FILE_WRITE_ERROR; 00354 sense[j] = ( bound_id < 0 ) ? -1 : 1; 00355 bound_ents[j] = new_faces[abs( bound_id ) - 1]; 00356 } 00357 00358 EntityType etype; 00359 result = readMeshIface->get_ordered_vertices( &bound_ents[0], sense, num_bounding, 3, &bound_verts[0], etype ); 00360 CHECK( "Failed in get_ordered_vertices." ); 00361 00362 // Make the element 00363 result = mdbImpl->create_element( etype, &bound_verts[0], CN::VerticesPerEntity( etype ), new_handle ); 00364 CHECK( "Failed to create region." ); 00365 00366 result = mdbImpl->add_entities( this_gent, &new_handle, 1 ); 00367 CHECK( "Failed to add region to geom set." ); 00368 00369 if( file_id_tag ) regions[i] = new_handle; 00370 00371 n = fscanf( file_ptr, "%d ", &dum_int ); 00372 CHECKN( 1 ); 00373 } // End of reading regions 00374 00375 if( file_id_tag ) 00376 { 00377 result = readMeshIface->assign_ids( *file_id_tag, ®ions[0], regions.size(), 1 ); 00378 if( MB_SUCCESS != result ) return result; 00379 } 00380 00381 return MB_SUCCESS; 00382 } 00383 00384 ErrorCode ReadSms::get_set( std::vector< EntityHandle >* sets, int set_dim, int set_id, Tag dim_tag, 00385 EntityHandle& this_set, const Tag* file_id_tag ) 00386 { 00387 ErrorCode result = MB_SUCCESS; 00388 00389 if( set_dim < 0 || set_dim > 3 ) return MB_FILE_WRITE_ERROR; 00390 00391 if( (int)sets[set_dim].size() <= set_id || !sets[set_dim][set_id] ) 00392 { 00393 if( (int)sets[set_dim].size() <= set_id ) sets[set_dim].resize( set_id + 1, 0 ); 00394 00395 if( !sets[set_dim][set_id] ) 00396 { 00397 result = mdbImpl->create_meshset( MESHSET_SET, sets[set_dim][set_id] ); 00398 if( MB_SUCCESS != result ) return result; 00399 result = mdbImpl->tag_set_data( globalId, &sets[set_dim][set_id], 1, &set_id ); 00400 if( MB_SUCCESS != result ) return result; 00401 result = mdbImpl->tag_set_data( dim_tag, &sets[set_dim][set_id], 1, &set_dim ); 00402 if( MB_SUCCESS != result ) return result; 00403 00404 if( file_id_tag ) 00405 { 00406 result = mdbImpl->tag_set_data( *file_id_tag, &sets[set_dim][set_id], 1, &setId ); 00407 ++setId; 00408 } 00409 } 00410 } 00411 00412 this_set = sets[set_dim][set_id]; 00413 00414 return result; 00415 } 00416 00417 ErrorCode ReadSms::read_parallel_info( FILE* file_ptr ) 00418 { 00419 // ErrorCode result; 00420 00421 // Read partition info 00422 int nparts, part_id, num_ifaces, num_corner_ents; 00423 int num_read = fscanf( file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents ); 00424 if( !num_read ) return MB_FAILURE; 00425 00426 // Read interfaces 00427 int iface_id, iface_dim, iface_own, num_iface_corners; 00428 // EntityHandle this_iface; 00429 std::vector< int >* iface_corners = NULL; 00430 for( int i = 0; i < num_ifaces; i++ ) 00431 { 00432 num_read = fscanf( file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own, &num_iface_corners ); 00433 if( !num_read ) return MB_FAILURE; 00434 00435 // result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface); 00436 // CHECK("Failed to make iface set."); 00437 00438 // Read the corner ids and store them on the set for now 00439 iface_corners = new std::vector< int >( num_iface_corners ); 00440 for( int j = 0; j < num_iface_corners; j++ ) 00441 { 00442 num_read = fscanf( file_ptr, "%d", &( *iface_corners )[j] ); 00443 if( !num_read ) 00444 { 00445 delete iface_corners; 00446 return MB_FAILURE; 00447 } 00448 } 00449 00450 // result = tag_set_data(ifaceCornerTag, &this_iface, 1, 00451 //&iface_corners); 00452 // CHECK("Failed to set iface corner tag."); 00453 00454 delete iface_corners; 00455 iface_corners = NULL; 00456 } 00457 00458 // Interface data has been read 00459 return MB_SUCCESS; 00460 } 00461 00462 ErrorCode ReadSms::add_entities( EntityHandle start, EntityHandle count, const Tag* file_id_tag ) 00463 { 00464 if( !count || !file_id_tag ) return MB_FAILURE; 00465 00466 Range range; 00467 range.insert( start, start + count - 1 ); 00468 return readMeshIface->assign_ids( *file_id_tag, range, 1 ); 00469 } 00470 00471 } // namespace moab