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