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