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