MOAB: Mesh Oriented datABase  (version 5.3.0)
ReadSms.cpp
Go to the documentation of this file.
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 <cerrno>
00033 #include <cstring>
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, &regions[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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines