MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 #include "ReadNASTRAN.hpp" 00017 00018 #include <iostream> 00019 #include <sstream> 00020 #include <fstream> 00021 #include <vector> 00022 #include <cstdlib> 00023 #include <assert.h> 00024 #include <cmath> 00025 00026 #include "moab/Interface.hpp" 00027 #include "moab/ReadUtilIface.hpp" 00028 #include "Internals.hpp" // For MB_START_ID 00029 #include "moab/Range.hpp" 00030 #include "moab/FileOptions.hpp" 00031 #include "FileTokenizer.hpp" 00032 #include "MBTagConventions.hpp" 00033 #include "moab/CN.hpp" 00034 00035 namespace moab 00036 { 00037 00038 ReaderIface* ReadNASTRAN::factory( Interface* iface ) 00039 { 00040 return new ReadNASTRAN( iface ); 00041 } 00042 00043 // Constructor 00044 ReadNASTRAN::ReadNASTRAN( Interface* impl ) : MBI( impl ) 00045 { 00046 assert( NULL != impl ); 00047 MBI->query_interface( readMeshIface ); 00048 assert( NULL != readMeshIface ); 00049 } 00050 00051 // Destructor 00052 ReadNASTRAN::~ReadNASTRAN() 00053 { 00054 if( readMeshIface ) 00055 { 00056 MBI->release_interface( readMeshIface ); 00057 readMeshIface = 0; 00058 } 00059 } 00060 00061 ErrorCode ReadNASTRAN::read_tag_values( const char* /*file_name*/, const char* /*tag_name*/, 00062 const FileOptions& /*opts*/, std::vector< int >& /*tag_values_out*/, 00063 const SubsetList* /*subset_list*/ ) 00064 { 00065 return MB_NOT_IMPLEMENTED; 00066 } 00067 00068 // Load the file as called by the Interface function 00069 ErrorCode ReadNASTRAN::load_file( const char* filename, const EntityHandle* /* file_set */, 00070 const FileOptions& /* opts */, const ReaderIface::SubsetList* subset_list, 00071 const Tag* file_id_tag ) 00072 { 00073 // At this time there is no support for reading a subset of the file 00074 if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for NASTRAN" ); } 00075 00076 nodeIdMap.clear(); 00077 elemIdMap.clear(); 00078 00079 bool debug = false; 00080 if( debug ) std::cout << "begin ReadNASTRAN::load_file" << std::endl; 00081 ErrorCode result; 00082 00083 // Count the entities of each type in the file. This is used to allocate the node array. 00084 int entity_count[MBMAXTYPE]; 00085 for( int i = 0; i < MBMAXTYPE; i++ ) 00086 entity_count[i] = 0; 00087 00088 /* Determine the line_format of the first line. Assume that the entire file 00089 has the same format. */ 00090 std::string line; 00091 std::ifstream file( filename ); 00092 if( !getline( file, line ) ) return MB_FILE_DOES_NOT_EXIST; 00093 line_format format; 00094 result = determine_line_format( line, format ); 00095 if( MB_SUCCESS != result ) return result; 00096 00097 /* Count the number of each entity in the file. This allows one to allocate 00098 a sequential array of vertex handles. */ 00099 while( !file.eof() ) 00100 { 00101 // Cut the line into fields as determined by the line format. 00102 // Use a vector to allow for an unknown number of tokens (continue lines). 00103 // Continue lines are not implemented. 00104 std::vector< std::string > tokens; 00105 tokens.reserve( 10 ); // assume 10 fields to avoid extra vector resizing 00106 result = tokenize_line( line, format, tokens ); 00107 if( MB_SUCCESS != result ) return result; 00108 00109 // Process the tokens of the line. The first token describes the entity type. 00110 EntityType type; 00111 result = determine_entity_type( ( tokens.empty() ) ? "" : tokens.front(), type ); 00112 if( MB_SUCCESS != result ) return result; 00113 entity_count[type]++; 00114 getline( file, line ); 00115 } 00116 00117 if( debug ) 00118 { 00119 for( int i = 0; i < MBMAXTYPE; i++ ) 00120 { 00121 std::cout << "entity_count[" << i << "]=" << entity_count[i] << std::endl; 00122 } 00123 } 00124 00125 // Keep list of material sets 00126 std::vector< Range > materials; 00127 00128 // Now that the number of vertices is known, create the vertices. 00129 EntityHandle start_vert = 0; 00130 std::vector< double* > coord_arrays( 3 ); 00131 result = readMeshIface->get_node_coords( 3, entity_count[0], MB_START_ID, start_vert, coord_arrays ); 00132 if( MB_SUCCESS != result ) return result; 00133 if( 0 == start_vert ) return MB_FAILURE; // check for NULL 00134 int id, vert_index = 0; 00135 if( debug ) std::cout << "allocated coord arrays" << std::endl; 00136 00137 // Read the file again to create the entities. 00138 file.clear(); // Clear eof state from object 00139 file.seekg( 0 ); // Rewind file 00140 while( !file.eof() ) 00141 { 00142 getline( file, line ); 00143 00144 // Cut the line into fields as determined by the line format. 00145 // Use a vector to allow for an unknown number of tokens (continue lines). 00146 // Continue lines are not implemented. 00147 std::vector< std::string > tokens; 00148 tokens.reserve( 10 ); // assume 10 fields to avoid extra vector resizing 00149 result = tokenize_line( line, format, tokens ); 00150 if( MB_SUCCESS != result ) return result; 00151 00152 // Process the tokens of the line. The first token describes the entity type. 00153 EntityType type; 00154 result = determine_entity_type( tokens.front(), type ); 00155 if( MB_SUCCESS != result ) return result; 00156 00157 // Create the entity. 00158 if( MBVERTEX == type ) 00159 { 00160 double* coords[3] = { coord_arrays[0] + vert_index, coord_arrays[1] + vert_index, 00161 coord_arrays[2] + vert_index }; 00162 result = read_node( tokens, debug, coords, id ); 00163 if( MB_SUCCESS != result ) return result; 00164 if( !nodeIdMap.insert( id, start_vert + vert_index, 1 ).second ) return MB_FAILURE; // Duplicate IDs! 00165 ++vert_index; 00166 } 00167 else 00168 { 00169 result = read_element( tokens, materials, type, debug ); 00170 if( MB_SUCCESS != result ) return result; 00171 } 00172 } 00173 00174 result = create_materials( materials ); 00175 if( MB_SUCCESS != result ) return result; 00176 00177 result = assign_ids( file_id_tag ); 00178 if( MB_SUCCESS != result ) return result; 00179 00180 file.close(); 00181 nodeIdMap.clear(); 00182 elemIdMap.clear(); 00183 return MB_SUCCESS; 00184 } 00185 00186 /* Determine the type of NASTRAN line: small field, large field, or free field. 00187 small field: each line has 10 fields of 8 characters 00188 large field: 1x8, 4x16, 1x8. Field 1 must have an asterisk following the character string 00189 free field: each line entry must be separated by a comma 00190 Implementation tries to avoid more searches than necessary. */ 00191 ErrorCode ReadNASTRAN::determine_line_format( const std::string& line, line_format& format ) 00192 { 00193 std::string::size_type found_asterisk = line.find( "*" ); 00194 if( std::string::npos != found_asterisk ) 00195 { 00196 format = LARGE_FIELD; 00197 return MB_SUCCESS; 00198 } 00199 else 00200 { 00201 std::string::size_type found_comma = line.find( "," ); 00202 if( std::string::npos != found_comma ) 00203 { 00204 format = FREE_FIELD; 00205 return MB_SUCCESS; 00206 } 00207 else 00208 { 00209 format = SMALL_FIELD; 00210 return MB_SUCCESS; 00211 } 00212 } 00213 } 00214 00215 /* Tokenize the line. Continue-lines have not been implemented. */ 00216 ErrorCode ReadNASTRAN::tokenize_line( const std::string& line, const line_format format, 00217 std::vector< std::string >& tokens ) 00218 { 00219 size_t line_size = line.size(); 00220 switch( format ) 00221 { 00222 case SMALL_FIELD: { 00223 // Expect 10 fields of 8 characters. 00224 // The sample file does not have all 10 fields in each line 00225 const int field_length = 8; 00226 unsigned int n_tokens = line_size / field_length; 00227 for( unsigned int i = 0; i < n_tokens; i++ ) 00228 { 00229 tokens.push_back( line.substr( i * field_length, field_length ) ); 00230 } 00231 break; 00232 } 00233 case LARGE_FIELD: 00234 return MB_NOT_IMPLEMENTED; 00235 case FREE_FIELD: 00236 return MB_NOT_IMPLEMENTED; 00237 default: 00238 return MB_FAILURE; 00239 } 00240 00241 return MB_SUCCESS; 00242 } 00243 00244 ErrorCode ReadNASTRAN::determine_entity_type( const std::string& first_token, EntityType& type ) 00245 { 00246 if( 0 == first_token.compare( "GRID " ) ) 00247 type = MBVERTEX; 00248 else if( 0 == first_token.compare( "CTETRA " ) ) 00249 type = MBTET; 00250 else if( 0 == first_token.compare( "CPENTA " ) ) 00251 type = MBPRISM; 00252 else if( 0 == first_token.compare( "CHEXA " ) ) 00253 type = MBHEX; 00254 else 00255 return MB_NOT_IMPLEMENTED; 00256 00257 return MB_SUCCESS; 00258 } 00259 00260 /* Some help from Jason: 00261 Nastran floats must contain a decimal point, may contain 00262 a leading '-' and may contain an exponent. The 'E' is optional 00263 when specifying an exponent. A '-' or '+' at any location other 00264 than the first position indicates an exponent. For a positive 00265 exponent, either a '+' or an 'E' must be specified. For a 00266 negative exponent, the 'E' is option and the '-' is always specified. 00267 Examples for the real value 7.0 from mcs2006 quick reference guide: 00268 7.0 .7E1 0.7+1 .70+1 7.E+0 70.-1 00269 00270 From the test file created in SC/Tetra: 00271 GRID 1 03.9804546.9052-15.6008-1 00272 has the coordinates: ( 3.980454, 6.9052e-1, 5.6008e-1 ) 00273 GRID 200005 04.004752-3.985-15.4955-1 00274 has the coordinates: ( 4.004752, -3.985e-1, 5.4955e-1 ) */ 00275 ErrorCode ReadNASTRAN::get_real( const std::string& token, double& real ) 00276 { 00277 std::string significand = token; 00278 std::string exponent = "0"; 00279 00280 // Cut off the first digit because a "-" could be here indicating a negative 00281 // number. Instead we are looking for a negative exponent. 00282 std::string back_token = token.substr( 1 ); 00283 00284 // A minus that is not the first digit is always a negative exponent 00285 std::string::size_type found_minus = back_token.find( "-" ); 00286 if( std::string::npos != found_minus ) 00287 { 00288 // separate the significand from the exponent at the "-" 00289 exponent = token.substr( found_minus + 1 ); 00290 significand = token.substr( 0, found_minus + 1 ); 00291 00292 // If the significand has an "E", remove it 00293 if( std::string::npos != significand.find( "E" ) ) 00294 // Assume the "E" is at the end of the significand. 00295 significand = significand.substr( 1, significand.size() - 2 ); 00296 00297 // If a minus does not exist past the 1st digit, but an "E" or "+" does, then 00298 // it is a positive exponent. First look for an "E", 00299 } 00300 else 00301 { 00302 std::string::size_type found_E = token.find( "E" ); 00303 if( std::string::npos != found_E ) 00304 { 00305 significand = token.substr( 0, found_E - 1 ); 00306 exponent = token.substr( found_E + 1 ); 00307 // If there is a "+" on the exponent, cut it off 00308 std::size_t found_plus = exponent.find( "+" ); 00309 if( std::string::npos != found_plus ) { exponent = exponent.substr( found_plus + 1 ); } 00310 } 00311 else 00312 { 00313 // If there is a "+" on the exponent, cut it off 00314 std::size_t found_plus = token.find( "+" ); 00315 if( std::string::npos != found_plus ) 00316 { 00317 significand = token.substr( 0, found_plus - 1 ); 00318 exponent = token.substr( found_plus + 1 ); 00319 } 00320 } 00321 } 00322 00323 // Now assemble the real number 00324 double signi = atof( significand.c_str() ); 00325 double expon = atof( exponent.c_str() ); 00326 00327 if( HUGE_VAL == signi || HUGE_VAL == expon ) return MB_FAILURE; 00328 00329 real = signi * pow( 10, expon ); 00330 00331 return MB_SUCCESS; 00332 } 00333 00334 /* It has been determined that this line is a vertex. Read the rest of 00335 the line and create the vertex. */ 00336 ErrorCode ReadNASTRAN::read_node( const std::vector< std::string >& tokens, const bool debug, double* coords[3], 00337 int& id ) 00338 { 00339 // Read the node's id (unique) 00340 ErrorCode result; 00341 id = atoi( tokens[1].c_str() ); 00342 00343 // Read the node's coordinate system number 00344 // "0" or blank refers to the basic coordinate system. 00345 int coord_system = atoi( tokens[2].c_str() ); 00346 if( 0 != coord_system ) 00347 { 00348 std::cerr << "ReadNASTRAN: alternative coordinate systems not implemented" << std::endl; 00349 return MB_NOT_IMPLEMENTED; 00350 } 00351 00352 // Read the coordinates 00353 for( unsigned int i = 0; i < 3; i++ ) 00354 { 00355 result = get_real( tokens[i + 3], *coords[i] ); 00356 if( MB_SUCCESS != result ) return result; 00357 if( debug ) std::cout << "read_node: coords[" << i << "]=" << coords[i] << std::endl; 00358 } 00359 00360 return MB_SUCCESS; 00361 } 00362 00363 /* The type of element has already been identified. Read the rest of the 00364 line and create the element. Assume that all of the nodes have already 00365 been created. */ 00366 ErrorCode ReadNASTRAN::read_element( const std::vector< std::string >& tokens, std::vector< Range >& materials, 00367 const EntityType element_type, const bool /*debug*/ ) 00368 { 00369 // Read the element's id (unique) and material set 00370 ErrorCode result; 00371 int id = atoi( tokens[1].c_str() ); 00372 int material = atoi( tokens[2].c_str() ); 00373 00374 // Resize materials list if necessary. This code is somewhat complicated 00375 // so as to avoid copying of Ranges 00376 if( material >= (int)materials.size() ) 00377 { 00378 if( (int)materials.capacity() < material ) 00379 materials.resize( material + 1 ); 00380 else 00381 { 00382 std::vector< Range > new_mat( material + 1 ); 00383 for( size_t i = 0; i < materials.size(); ++i ) 00384 new_mat[i].swap( materials[i] ); 00385 materials.swap( new_mat ); 00386 } 00387 } 00388 00389 // The size of the connectivity array depends on the element type 00390 int n_conn = CN::VerticesPerEntity( element_type ); 00391 EntityHandle conn_verts[27]; 00392 assert( n_conn <= (int)( sizeof( conn_verts ) / sizeof( EntityHandle ) ) ); 00393 00394 // Read the connected node ids from the file 00395 for( int i = 0; i < n_conn; i++ ) 00396 { 00397 int n = atoi( tokens[3 + i].c_str() ); 00398 conn_verts[i] = nodeIdMap.find( n ); 00399 if( !conn_verts[i] ) // invalid vertex id 00400 return MB_FAILURE; 00401 } 00402 00403 // Create the element and set the global id from the NASTRAN file 00404 EntityHandle element; 00405 result = MBI->create_element( element_type, conn_verts, n_conn, element ); 00406 if( MB_SUCCESS != result ) return result; 00407 elemIdMap.insert( id, element, 1 ); 00408 00409 materials[material].insert( element ); 00410 return MB_SUCCESS; 00411 } 00412 00413 ErrorCode ReadNASTRAN::create_materials( const std::vector< Range >& materials ) 00414 { 00415 ErrorCode result; 00416 Tag material_tag; 00417 int negone = -1; 00418 result = MBI->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, material_tag, MB_TAG_SPARSE | MB_TAG_CREAT, 00419 &negone ); 00420 if( MB_SUCCESS != result ) return result; 00421 00422 for( size_t i = 0; i < materials.size(); ++i ) 00423 { 00424 if( materials[i].empty() ) continue; 00425 00426 // Merge with existing or create new? Original code effectively 00427 // created new by only merging with existing in current file set, 00428 // so do the same here. - j.kraftcheck 00429 00430 EntityHandle handle; 00431 result = MBI->create_meshset( MESHSET_SET, handle ); 00432 if( MB_SUCCESS != result ) return result; 00433 00434 result = MBI->add_entities( handle, materials[i] ); 00435 if( MB_SUCCESS != result ) return result; 00436 00437 int id = i; 00438 result = MBI->tag_set_data( material_tag, &handle, 1, &id ); 00439 if( MB_SUCCESS != result ) return result; 00440 } 00441 00442 return MB_SUCCESS; 00443 } 00444 00445 ErrorCode ReadNASTRAN::assign_ids( const Tag* file_id_tag ) 00446 { 00447 // Create tag 00448 ErrorCode result; 00449 Tag id_tag = MBI->globalId_tag(); 00450 00451 RangeMap< int, EntityHandle >::iterator i; 00452 for( int t = 0; t < 2; ++t ) 00453 { 00454 RangeMap< int, EntityHandle >& fileIdMap = t ? elemIdMap : nodeIdMap; 00455 for( i = fileIdMap.begin(); i != fileIdMap.end(); ++i ) 00456 { 00457 Range range( i->value, i->value + i->count - 1 ); 00458 00459 result = readMeshIface->assign_ids( id_tag, range, i->begin ); 00460 if( MB_SUCCESS != result ) return result; 00461 00462 if( file_id_tag && *file_id_tag != id_tag ) 00463 { 00464 result = readMeshIface->assign_ids( *file_id_tag, range, i->begin ); 00465 if( MB_SUCCESS != result ) return result; 00466 } 00467 } 00468 } 00469 00470 return MB_SUCCESS; 00471 } 00472 00473 } // namespace moab