MOAB: Mesh Oriented datABase  (version 5.3.1)
ReadTetGen.cpp
Go to the documentation of this file.
00001 #include "ReadTetGen.hpp"
00002 #include "moab/Interface.hpp"
00003 #include "moab/Range.hpp"
00004 #include "moab/ReadUtilIface.hpp"
00005 #include "moab/FileOptions.hpp"
00006 #include "MBTagConventions.hpp"
00007 #include <iostream>
00008 #include <fstream>
00009 #include <sstream>
00010 #include <cctype>
00011 #include <map>
00012 
00013 namespace moab
00014 {
00015 
00016 ReaderIface* ReadTetGen::factory( Interface* moab )
00017 {
00018     return new ReadTetGen( moab );
00019 }
00020 
00021 ReadTetGen::ReadTetGen( Interface* moab ) : mbIface( moab ), readTool( 0 )
00022 {
00023     moab->query_interface( readTool );
00024 }
00025 
00026 ReadTetGen::~ReadTetGen()
00027 {
00028     if( mbIface && readTool ) mbIface->release_interface( readTool );
00029 }
00030 
00031 ErrorCode ReadTetGen::open_file( const std::string& filename, const std::string& basename, const std::string& suffix,
00032                                  const char* exp_suffix, const char* opt_name, const FileOptions& opts,
00033                                  std::ifstream& file_stream, bool file_required )
00034 {
00035     std::string real_file_name;
00036     ErrorCode rval = opts.get_option( opt_name, real_file_name );
00037     if( MB_ENTITY_NOT_FOUND == rval || real_file_name.empty() )
00038     {
00039         if( MB_SUCCESS == rval ) file_required = true;
00040         if( suffix == exp_suffix ) { real_file_name = filename; }
00041         else
00042         {
00043             real_file_name = basename;
00044             real_file_name += ".";
00045             real_file_name += exp_suffix;
00046         }
00047     }
00048 
00049     if( !real_file_name.empty() ) file_stream.open( real_file_name.c_str(), std::ios::in );
00050     if( file_required && !file_stream.is_open() )
00051     {
00052         MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, real_file_name << ": cannot read file" );
00053     }
00054 
00055     return MB_SUCCESS;
00056 }
00057 
00058 ErrorCode ReadTetGen::read_tag_values( const char* /* file_name */, const char* /* tag_name */,
00059                                        const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */,
00060                                        const SubsetList* /* subset_list */ )
00061 {
00062     return MB_NOT_IMPLEMENTED;
00063 }
00064 
00065 ErrorCode ReadTetGen::load_file( const char* file_name_c, const EntityHandle* /* file_set */, const FileOptions& opts,
00066                                  const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag )
00067 {
00068     std::ifstream node_file, ele_file, face_file, edge_file;
00069     ErrorCode rval;
00070 
00071     if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for TetGen" ); }
00072 
00073     std::string suffix, base, filename( file_name_c );
00074     size_t dot_idx = filename.find_last_of( '.' );
00075     if( dot_idx == std::string::npos ) { base = filename; }
00076     else
00077     {
00078         suffix = filename.substr( dot_idx + 1 );
00079         for( size_t i = 0; i < suffix.length(); ++i )
00080             suffix[i] = (char)tolower( suffix[i] );
00081         if( suffix == "node" || suffix == "ele" || suffix == "face" || suffix == "edge" )
00082         {
00083             base = filename.substr( 0, dot_idx );
00084         }
00085         else
00086         {
00087             base = filename;
00088             suffix.clear();
00089         }
00090     }
00091 
00092     rval = open_file( filename, base, suffix, "node", "NODE_FILE", opts, node_file, true );
00093     if( MB_SUCCESS != rval ) return rval;
00094     rval = open_file( filename, base, suffix, "ele", "ELE_FILE", opts, ele_file );
00095     if( MB_SUCCESS != rval ) return rval;
00096     rval = open_file( filename, base, suffix, "face", "FACE_FILE", opts, face_file );
00097     if( MB_SUCCESS != rval ) return rval;
00098     rval = open_file( filename, base, suffix, "edge", "EDGE_FILE", opts, edge_file );
00099     if( MB_SUCCESS != rval ) return rval;
00100 
00101     std::vector< Tag > attr_tags[4];
00102     std::vector< int > attr_idx[4];
00103     const char* option_names[4] = { "NODE_ATTR_LIST", "EDGE_ATTR_LIST", "TRI_ATTR_LIST", "TET_ATTR_LIST" };
00104     const char* group_names[4]  = { 0, "CURVE_ID", "SURFACE_ID", "VOLUME_ID" };
00105     for( int i = 0; i < 4; ++i )
00106     {
00107         std::string opt_str;
00108         rval = opts.get_str_option( option_names[i], opt_str );
00109         if( MB_SUCCESS != rval ) continue;
00110         rval = parse_attr_list( opt_str, attr_tags[i], attr_idx[i], group_names[i] );
00111         if( MB_SUCCESS != rval ) { MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, option_names[i] << ": invalid option value" ); }
00112     }
00113 
00114     Range tets, tris, edges;
00115     std::vector< EntityHandle > nodes;
00116     rval = read_node_file( node_file, &attr_tags[0][0], &attr_idx[0][0], attr_tags[0].size(), nodes );
00117     if( MB_SUCCESS == rval && ele_file.is_open() ) rval = read_elem_file( MBTET, ele_file, nodes, tets );
00118     if( MB_SUCCESS == rval && face_file.is_open() ) rval = read_elem_file( MBTRI, face_file, nodes, tris );
00119     if( MB_SUCCESS == rval && edge_file.is_open() ) rval = read_elem_file( MBEDGE, edge_file, nodes, edges );
00120 
00121     if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, &nodes[0], nodes.size() );
00122     if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, edges );
00123     if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tris );
00124     if( file_id_tag && MB_SUCCESS == rval ) rval = readTool->assign_ids( *file_id_tag, tets );
00125 
00126     return rval;
00127 }
00128 
00129 ErrorCode ReadTetGen::parse_attr_list( const std::string& option_str, std::vector< Tag >& tag_list,
00130                                        std::vector< int >& index_list, const char* group_designator )
00131 {
00132     std::vector< std::string > name_list;
00133     size_t prev_pos = 0;
00134     while( prev_pos != std::string::npos )
00135     {
00136         size_t pos = option_str.find_first_of( ',', prev_pos );
00137         name_list.push_back( option_str.substr( prev_pos, pos ) );
00138         prev_pos = pos + 1;
00139     }
00140 
00141     index_list.resize( name_list.size() );
00142     std::map< std::string, int > name_count;
00143     for( size_t i = 0; i < name_list.size(); ++i )
00144         index_list[i] = name_count[name_list[i]]++;
00145 
00146     for( size_t i = 0; i < name_list.size(); ++i )
00147     {
00148         if( group_designator && name_list[i] == group_designator )
00149         {
00150             tag_list[i]   = 0;
00151             index_list[i] = -1;
00152         }
00153         else if( name_list.empty() )
00154         {
00155             tag_list[i]   = 0;
00156             index_list[i] = 0;
00157         }
00158         else
00159         {
00160             ErrorCode rval = mbIface->tag_get_handle( name_list[i].c_str(), name_count[name_list[i]], MB_TYPE_DOUBLE,
00161                                                       tag_list[i], MB_TAG_DENSE | MB_TAG_CREAT );
00162             if( MB_SUCCESS != rval ) return rval;
00163         }
00164     }
00165 
00166     return MB_SUCCESS;
00167 }
00168 
00169 ErrorCode ReadTetGen::read_line( std::istream& file, std::string& line, int& lineno )
00170 {
00171     // Loop until we find a non-empty line
00172     do
00173     {
00174         // Read a line
00175         line.clear();
00176         if( !getline( file, line ) ) return MB_FILE_WRITE_ERROR;
00177         ++lineno;
00178         // Strip comments from line
00179         size_t pos = line.find_first_of( '#' );
00180         if( pos != std::string::npos ) line = line.substr( 0, pos );
00181         // Strip leading whitespace from line
00182         for( pos = 0; pos < line.length() && isspace( line[pos] ); ++pos )
00183             ;
00184         if( pos == line.length() )
00185             line.clear();
00186         else if( pos != 0 )
00187             line = line.substr( pos );
00188     } while( line.empty() );
00189 
00190     return MB_SUCCESS;
00191 }
00192 
00193 ErrorCode ReadTetGen::read_line( std::istream& file, double* values_out, int num_values, int& lineno )
00194 {
00195     // Get a line of text
00196     std::string line;
00197     ErrorCode rval = read_line( file, line, lineno );
00198     if( MB_SUCCESS != rval ) return rval;
00199 
00200     // Tokenize line as doubles
00201     std::stringstream str( line );
00202     for( int i = 0; i < num_values; ++i )
00203     {
00204         double v;
00205         if( !( str >> v ) ) { MB_SET_ERR( MB_FAILURE, "Error reading node data at line " << lineno ); }
00206         values_out[i] = v;
00207     }
00208 
00209     // Check that we're at the end of the line
00210     int junk;
00211     if( ( str >> junk ) || !str.eof() )
00212     {
00213         MB_SET_ERR( MB_FAILURE, "Unexpected trailing data for line " << lineno << " of node data" );
00214     }
00215 
00216     return MB_SUCCESS;
00217 }
00218 
00219 ErrorCode ReadTetGen::read_node_file( std::istream& file, const Tag* attr_tag_list, const int* attr_tag_index,
00220                                       int attr_tag_list_len, std::vector< EntityHandle >& nodes )
00221 {
00222     int lineno = 0;
00223     ErrorCode rval;
00224 
00225     double header_vals[4];
00226     rval = read_line( file, header_vals, 4, lineno );
00227     if( MB_SUCCESS != rval ) return rval;
00228 
00229     const int num_vtx   = (int)header_vals[0];
00230     const int dim       = (int)header_vals[1];
00231     const int num_attr  = (int)header_vals[2];
00232     const int bdry_flag = (int)header_vals[3];
00233     if( num_vtx < 1 || dim < 2 || dim > 3 || num_attr < 0 || bdry_flag < 0 || bdry_flag > 1 )
00234     {
00235         MB_SET_ERR( MB_FAILURE, "Invalid header line for node data" );
00236     }
00237     if( attr_tag_list_len > num_attr ) attr_tag_list_len = num_attr;
00238 
00239     // Allocate space for tag data
00240     std::map< Tag, int > tag_size;
00241     std::map< Tag, std::vector< double > > tag_data;
00242     for( int i = 0; i < attr_tag_list_len; ++i )
00243     {
00244         if( !attr_tag_list[i] || attr_tag_index[i] < 0 ) continue;
00245         std::vector< double >& data = tag_data[attr_tag_list[i]];
00246         // Increase tag size by one value per vertex for each time
00247         // we encounter it in the list.
00248         data.resize( data.size() + num_vtx );
00249         ++tag_size[attr_tag_list[i]];
00250     }
00251     std::vector< double* > attr_data( attr_tag_list_len );
00252     std::vector< int > attr_size( attr_tag_list_len );
00253     for( int i = 0; i < attr_tag_list_len; ++i )
00254     {
00255         if( !attr_tag_list[i] || attr_tag_index[i] < 0 )
00256         {
00257             attr_data[i] = 0;
00258             attr_size[i] = 0;
00259         }
00260         else
00261         {
00262             attr_data[i] = &( tag_data[attr_tag_list[i]] )[0];
00263             attr_size[i] = tag_size[attr_tag_list[i]];
00264         }
00265     }
00266 
00267     // Allocate vertices
00268     std::vector< double* > coords;
00269     EntityHandle start_handle;
00270     rval = readTool->get_node_coords( dim, num_vtx, 1, start_handle, coords );
00271     if( MB_SUCCESS != rval ) return rval;
00272 
00273     // Read data line for each node
00274     nodes.reserve( num_vtx );
00275     std::vector< double > data( 1 + dim + num_attr + bdry_flag );
00276     std::vector< int > ids( num_vtx );
00277     for( int i = 0; i < num_vtx; ++i )
00278     {
00279         rval = read_line( file, &data[0], data.size(), lineno );
00280         if( MB_SUCCESS != rval ) return rval;
00281 
00282         // Get ID
00283         ids[i] = (int)data[0];
00284         if( ids[i] >= (int)nodes.size() ) nodes.resize( ids[i] + 1 );
00285         nodes[ids[i]] = start_handle + i;
00286 
00287         // Get coordinates
00288         // Cppcheck warning (false positive): variable coords is assigned a value that is never used
00289         for( int j = 0; j < dim; ++j )
00290             coords[j][i] = data[j + 1];
00291 
00292         // Get attribute data
00293         for( int j = 0; j < attr_tag_list_len; ++j )
00294             if( attr_data[j] ) attr_data[j][i * attr_size[j] + attr_tag_index[j]] = data[j + 1 + dim];
00295 
00296         // Discard boundary bit
00297     }
00298 
00299     // Store tag data
00300     Range node_range;
00301     node_range.insert( start_handle, start_handle + num_vtx - 1 );
00302     for( std::map< Tag, std::vector< double > >::iterator i = tag_data.begin(); i != tag_data.end(); ++i )
00303     {
00304         rval = mbIface->tag_set_data( i->first, node_range, &i->second[0] );
00305         if( MB_SUCCESS != rval ) return rval;
00306     }
00307     Tag idtag = mbIface->globalId_tag();
00308 
00309     rval = mbIface->tag_set_data( idtag, node_range, &ids[0] );
00310     if( MB_SUCCESS != rval ) return rval;
00311 
00312     return MB_SUCCESS;
00313 }
00314 
00315 ErrorCode ReadTetGen::read_elem_file( EntityType type, std::istream& file, const std::vector< EntityHandle >& nodes,
00316                                       Range& elems )
00317 {
00318     int lineno = 0;
00319     ErrorCode rval;
00320 
00321     int node_per_elem, have_group_id, dim;
00322     double header_vals[3];
00323     switch( type )
00324     {
00325         case MBTET:
00326             rval          = read_line( file, header_vals, 3, lineno );
00327             node_per_elem = (int)header_vals[1];
00328             have_group_id = (int)header_vals[2];
00329             dim           = 3;
00330             break;
00331         case MBTRI:
00332             rval          = read_line( file, header_vals, 2, lineno );
00333             node_per_elem = 3;
00334             have_group_id = (int)header_vals[1];
00335             dim           = 2;
00336             break;
00337         case MBEDGE:
00338             rval          = read_line( file, header_vals, 1, lineno );
00339             node_per_elem = 2;
00340             have_group_id = 0;
00341             dim           = 1;
00342             break;
00343         default:
00344             rval = MB_FAILURE;
00345             break;
00346     }
00347     if( MB_SUCCESS != rval ) return rval;
00348     const int num_elem = (int)header_vals[0];
00349     if( num_elem < 1 || node_per_elem < 2 || have_group_id < 0 || have_group_id > 1 )
00350     {
00351         MB_SET_ERR( MB_FAILURE, "Invalid header line for element data" );
00352     }
00353 
00354     // Create group map
00355     std::map< double, EntityHandle > groups;
00356     Tag dim_tag, id_tag;
00357     id_tag = mbIface->globalId_tag();
00358 
00359     const int negone = -1;
00360     rval = mbIface->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
00361                                     &negone );
00362     if( MB_SUCCESS != rval ) return rval;
00363 
00364     // Allocate elements
00365     EntityHandle start_handle, *conn_array;
00366     rval = readTool->get_element_connect( num_elem, node_per_elem, type, 1, start_handle, conn_array );
00367     if( MB_SUCCESS != rval ) return rval;
00368     elems.insert( start_handle, start_handle + num_elem - 1 );
00369 
00370     // Read data line for each node
00371     std::vector< double > data( 1 + node_per_elem + have_group_id );
00372     std::vector< int > ids( num_elem );
00373     for( int i = 0; i < num_elem; ++i )
00374     {
00375         rval = read_line( file, &data[0], data.size(), lineno );
00376         if( MB_SUCCESS != rval ) return rval;
00377 
00378         // Get ID
00379         ids[i] = (int)data[0];
00380 
00381         // Get connectivity
00382         for( int j = 0; j < node_per_elem; ++j )
00383             conn_array[node_per_elem * i + j] = nodes[(int)data[j + 1]];
00384 
00385         // Grouping
00386         if( have_group_id && 0.0 != data[node_per_elem + 1] )
00387         {
00388             double id        = data[node_per_elem + 1];
00389             EntityHandle grp = groups[id];
00390             if( 0 == grp )
00391             {
00392                 rval = mbIface->create_meshset( MESHSET_SET, grp );
00393                 if( MB_SUCCESS != rval ) return rval;
00394                 elems.insert( grp );
00395                 rval = mbIface->tag_set_data( dim_tag, &grp, 1, &dim );
00396                 if( MB_SUCCESS != rval ) return rval;
00397                 int iid = (int)id;
00398                 rval    = mbIface->tag_set_data( id_tag, &grp, 1, &iid );
00399                 if( MB_SUCCESS != rval ) return rval;
00400                 groups[id] = grp;
00401             }
00402             EntityHandle handle = start_handle + i;
00403             rval                = mbIface->add_entities( grp, &handle, 1 );
00404             if( MB_SUCCESS != rval ) return rval;
00405         }
00406     }
00407 
00408     // Store id data
00409     Range elems2;
00410     elems2.insert( start_handle, start_handle + num_elem - 1 );
00411     rval = mbIface->tag_set_data( id_tag, elems2, &ids[0] );
00412     if( MB_SUCCESS != rval ) return rval;
00413 
00414     return MB_SUCCESS;
00415 }
00416 
00417 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines