MOAB: Mesh Oriented datABase  (version 5.4.0)
ReadIDEAS.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <vector>
00004 #include <cstdlib>
00005 #include <sstream>
00006 #include <cassert>
00007 
00008 #include "ReadIDEAS.hpp"
00009 #include "moab/Interface.hpp"
00010 #include "Internals.hpp"
00011 #include "moab/ReadUtilIface.hpp"
00012 #include "FileTokenizer.hpp"
00013 #include "MBTagConventions.hpp"
00014 #include "moab/Range.hpp"
00015 #include "moab/CN.hpp"
00016 
00017 namespace moab
00018 {
00019 
00020 ReaderIface* ReadIDEAS::factory( Interface* iface )
00021 {
00022     return new ReadIDEAS( iface );
00023 }
00024 
00025 ReadIDEAS::ReadIDEAS( Interface* impl ) : MBI( impl )
00026 {
00027     impl->query_interface( readMeshIface );
00028 }
00029 
00030 ErrorCode ReadIDEAS::read_tag_values( const char* /* file_name */,
00031                                       const char* /* tag_name */,
00032                                       const FileOptions& /* opts */,
00033                                       std::vector< int >& /* tag_values_out */,
00034                                       const SubsetList* /* subset_list */ )
00035 {
00036     return MB_NOT_IMPLEMENTED;
00037 }
00038 
00039 ErrorCode ReadIDEAS::load_file( const char* fname,
00040                                 const EntityHandle*,
00041                                 const FileOptions& /*options*/,
00042                                 const ReaderIface::SubsetList* subset_list,
00043                                 const Tag* file_id_tag )
00044 {
00045     if( subset_list )
00046     {
00047         MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for IDEAS" );
00048     }
00049 
00050     file.open( fname );
00051     if( !file.good() )
00052     {
00053         MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "Failed to open file: " << fname );
00054     }
00055 
00056     ErrorCode rval;
00057 
00058     char line[10000];
00059     file.getline( line, 10000 );
00060     char* liter = line;
00061     while( *liter && isspace( *liter ) )
00062         ++liter;
00063     if( *liter != '-' ) return MB_FAILURE;
00064     ++liter;
00065     if( *liter != '1' ) return MB_FAILURE;
00066     while( *++liter )
00067         if( !isspace( *liter ) ) return MB_FAILURE;
00068 
00069     EntityHandle first_vertex = 0;
00070     while( !file.eof() )
00071     {
00072         file.getline( line, 10000 );
00073         unsigned int header_id = (unsigned int)strtol( line, NULL, 10 );
00074 
00075         // Create vertices
00076         if( DOUBLE_PRECISION_NODES0 == header_id || DOUBLE_PRECISION_NODES1 == header_id )
00077         {
00078             if( first_vertex )  // multiple vertex blocks?
00079                 return MB_FAILURE;
00080             rval = create_vertices( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read vertices" );
00081         }
00082         // Create elements
00083         else if( ELEMENTS0 == header_id || ELEMENTS1 == header_id || ELEMENTS2 == header_id )
00084         {
00085             if( !first_vertex )  // Need to read vertices first
00086                 return MB_FAILURE;
00087             rval = create_elements( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read elements" );
00088         }
00089         // Skip everything else
00090         else
00091         {
00092             rval = skip_header();
00093             if( MB_SUCCESS != rval ) return MB_FAILURE;
00094         }
00095     }
00096 
00097     file.close();
00098     return MB_SUCCESS;
00099 }
00100 
00101 ErrorCode ReadIDEAS::skip_header()
00102 {
00103     // Go until finding a pair of -1 lines
00104     char* ctmp;
00105     char line[10000];
00106     std::string s;
00107 
00108     int end_of_block = 0;
00109 
00110     long int il;
00111 
00112     while( file.getline( line, 10000 ) )
00113     {
00114         il = std::strtol( line, &ctmp, 10 );
00115         if( il == -1 )
00116         {
00117             s = ctmp;
00118             if( s.empty() ) end_of_block++;
00119         }
00120         else
00121             end_of_block = 0;
00122 
00123         if( end_of_block >= 2 ) return MB_SUCCESS;
00124     }
00125 
00126     return MB_SUCCESS;
00127 }
00128 
00129 ErrorCode ReadIDEAS::create_vertices( EntityHandle& first_vertex, const Tag* file_id_tag )
00130 {
00131     // Read two lines: first has some data, second has coordinates
00132     char line1[10000], line2[10000];
00133     int il1, il2;
00134     char *ctmp1, *ctmp2;
00135     std::string s1, s2;
00136 
00137     ErrorCode rval;
00138 
00139     std::streampos top_of_block = file.tellg();
00140     unsigned int num_verts      = 0;
00141 
00142     for( ;; )
00143     {
00144         // Read both lines
00145         if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
00146         if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
00147 
00148         // Check if we are at the end of the block
00149         il1 = std::strtol( line1, &ctmp1, 10 );
00150         il2 = std::strtol( line2, &ctmp2, 10 );
00151         if( ( il1 == -1 ) && ( il2 == -1 ) )
00152         {
00153             s1 = ctmp1;
00154             s2 = ctmp2;
00155             if( ( s1.empty() ) && ( s2.empty() ) ) break;
00156         }
00157         num_verts++;
00158     }
00159 
00160     file.seekg( top_of_block );
00161 
00162     std::vector< double* > arrays;
00163     rval = readMeshIface->get_node_coords( 3, num_verts, 0, first_vertex, arrays );
00164     if( MB_SUCCESS != rval ) return rval;
00165 
00166     Range verts;
00167     verts.insert( first_vertex, first_vertex + num_verts - 1 );
00168 
00169     double* x = arrays[0];
00170     double* y = arrays[1];
00171     double* z = arrays[2];
00172 
00173     // For now, assume ids are sequential and begin with 1
00174     Tag id_tag                  = MBI->globalId_tag();
00175     const int beginning_node_id = 1;
00176     int node_id                 = beginning_node_id;
00177 
00178     for( unsigned int i = 0; i < num_verts; i++ )
00179     {
00180         if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
00181         if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
00182 
00183         // Get the id out of the 1st line. Check the assumption that node ids are
00184         // sequential and begin with 1.
00185         if( node_id != std::strtol( line1, &ctmp1, 10 ) )
00186             MB_SET_ERR( MB_FAILURE, "node_id " << node_id << " line2:" << line2 << " ctmp1:" << ctmp1 );
00187         else
00188             ++node_id;
00189 
00190         // Get the doubles out of the 2nd line
00191         x[i] = std::strtod( line2, &ctmp2 );
00192         y[i] = std::strtod( ctmp2 + 1, &ctmp2 );
00193         z[i] = std::strtod( ctmp2 + 1, NULL );
00194     }
00195 
00196     if( !file.getline( line1, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines" );
00197     if( !file.getline( line2, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines 2" );
00198 
00199     // Tag the nodes with ids
00200     rval = readMeshIface->assign_ids( id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
00201     if( file_id_tag )
00202     {
00203         rval = readMeshIface->assign_ids( *file_id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
00204     }
00205 
00206     return MB_SUCCESS;
00207 }
00208 
00209 ErrorCode ReadIDEAS::create_elements( EntityHandle vstart, const Tag* file_id_tag )
00210 {
00211     char line1[10000], line2[10000];
00212     int il1, il2;
00213     char *ctmp1, *ctmp2;
00214     std::string s1, s2;
00215     ErrorCode rval;
00216     EntityHandle handle;
00217 
00218     Tag mat_tag, phys_tag, id_tag;
00219     rval = MBI->tag_get_handle( MAT_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, mat_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00220     if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
00221     rval = MBI->tag_get_handle( PHYS_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, phys_tag, MB_TAG_DENSE | MB_TAG_CREAT );
00222     if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
00223     id_tag = MBI->globalId_tag();
00224 
00225     for( ;; )
00226     {
00227         if( !file.getline( line1, 10000 ) || !file.getline( line2, 10000 ) ) return MB_FAILURE;
00228 
00229         // Check if we are at the end of the block
00230         il1 = std::strtol( line1, &ctmp1, 10 );
00231         il2 = std::strtol( line2, &ctmp2, 10 );
00232         if( ( il1 == -1 ) && ( il2 == -1 ) )
00233         {
00234             s1 = ctmp1;
00235             s2 = ctmp2;
00236             if( ( s1.empty() ) && ( s2.empty() ) ) return MB_SUCCESS;
00237         }
00238 
00239         // The first line describes attributes of the element other than connectivity.
00240         const int element_id = strtol( line1 + 1, &ctmp1, 10 );
00241         const int ideas_type = strtol( line1 + 11, &ctmp1, 10 );
00242         const int phys_table = strtol( line1 + 21, &ctmp1, 10 );
00243         const int mat_table  = strtol( line1 + 31, &ctmp1, 10 );
00244 
00245         // Determine the element type.
00246         EntityType mb_type;
00247         if( TRI0 == ideas_type || TRI1 == ideas_type )
00248             mb_type = MBTRI;
00249         else if( QUAD0 == ideas_type || QUAD1 == ideas_type )
00250             mb_type = MBQUAD;
00251         else if( TET == ideas_type )
00252             mb_type = MBTET;
00253         else if( HEX == ideas_type )
00254             mb_type = MBHEX;
00255         else if( WEDGE == ideas_type )
00256             mb_type = MBPRISM;
00257         else
00258         {
00259             std::cout << "IDEAS element type not yet added to MOAB reader." << std::endl;
00260             return MB_NOT_IMPLEMENTED;
00261         }
00262 
00263         // Get the connectivity out of the 2nd line
00264         std::stringstream ss( line2 );
00265         const int n_conn = CN::VerticesPerEntity( mb_type );
00266         EntityHandle conn[CN::MAX_NODES_PER_ELEMENT];
00267         EntityHandle vert;
00268         for( int i = 0; i < n_conn; ++i )
00269         {
00270             ss >> vert;
00271             conn[i] = vstart + vert - 1;
00272         }
00273 
00274         // Make the element. According to the Gmsh 2.2.3 source code, the IDEAS
00275         // canonical numbering is the same as MBCN.
00276         rval = MBI->create_element( mb_type, conn, n_conn, handle );MB_CHK_SET_ERR( rval, "can't create elements of type " << mb_type );
00277 
00278         // If the phys set does not already exist, create it.
00279         Range phys_sets;
00280         EntityHandle phys_set;
00281         const void* const phys_set_id_val[] = { &phys_table };
00282         rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &phys_tag, phys_set_id_val, 1, phys_sets );MB_CHK_SET_ERR( rval, "can't get phys sets" );
00283         if( phys_sets.empty() )
00284         {
00285             rval = MBI->create_meshset( MESHSET_SET, phys_set );MB_CHK_SET_ERR( rval, "can't create phys set" );
00286             rval = MBI->tag_set_data( phys_tag, &phys_set, 1, &phys_table );MB_CHK_SET_ERR( rval, "can't set tag to phys set" );
00287         }
00288         else if( 1 == phys_sets.size() )
00289         {
00290             phys_set = phys_sets.front();
00291         }
00292         else
00293         {
00294             return MB_MULTIPLE_ENTITIES_FOUND;
00295         }
00296         rval = MBI->add_entities( phys_set, &handle, 1 );MB_CHK_SET_ERR( rval, "can't add entities to phys set" );
00297 
00298         // If the material set does not already exist, create it.
00299         Range mat_sets;
00300         EntityHandle mat_set;
00301         const void* const mat_set_id_val[] = { &mat_table };
00302         rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, mat_set_id_val, 1, mat_sets );
00303         if( MB_SUCCESS != rval ) return rval;
00304         if( mat_sets.empty() )
00305         {
00306             rval = MBI->create_meshset( MESHSET_SET, mat_set );
00307             if( MB_SUCCESS != rval ) return rval;
00308             rval = MBI->tag_set_data( mat_tag, &mat_set, 1, &mat_table );
00309             if( MB_SUCCESS != rval ) return rval;
00310         }
00311         else if( 1 == mat_sets.size() )
00312         {
00313             mat_set = mat_sets.front();
00314         }
00315         else
00316         {
00317             return MB_MULTIPLE_ENTITIES_FOUND;
00318         }
00319         rval = MBI->add_entities( mat_set, &handle, 1 );
00320         if( MB_SUCCESS != rval ) return rval;
00321 
00322         // Tag the element with its id
00323         rval = MBI->tag_set_data( id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
00324         if( file_id_tag )
00325         {
00326             rval = MBI->tag_set_data( *file_id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
00327         }
00328     }
00329 
00330     return MB_SUCCESS;
00331 }
00332 
00333 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines