Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadNASTRAN.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #include "ReadNASTRAN.hpp"
17 
18 #include <iostream>
19 #include <sstream>
20 #include <fstream>
21 #include <vector>
22 #include <cstdlib>
23 #include <cassert>
24 #include <cmath>
25 
26 #include "moab/Interface.hpp"
27 #include "moab/ReadUtilIface.hpp"
28 #include "Internals.hpp" // For MB_START_ID
29 #include "moab/Range.hpp"
30 #include "moab/FileOptions.hpp"
31 #include "FileTokenizer.hpp"
32 #include "MBTagConventions.hpp"
33 #include "moab/CN.hpp"
34 
35 namespace moab
36 {
37 
39 {
40  return new ReadNASTRAN( iface );
41 }
42 
43 // Constructor
45 {
46  assert( NULL != impl );
48  assert( NULL != readMeshIface );
49 }
50 
51 // Destructor
53 {
54  if( readMeshIface )
55  {
57  readMeshIface = 0;
58  }
59 }
60 
61 ErrorCode ReadNASTRAN::read_tag_values( const char* /*file_name*/,
62  const char* /*tag_name*/,
63  const FileOptions& /*opts*/,
64  std::vector< int >& /*tag_values_out*/,
65  const SubsetList* /*subset_list*/ )
66 {
67  return MB_NOT_IMPLEMENTED;
68 }
69 
70 // Load the file as called by the Interface function
71 ErrorCode ReadNASTRAN::load_file( const char* filename,
72  const EntityHandle* /* file_set */,
73  const FileOptions& /* opts */,
74  const ReaderIface::SubsetList* subset_list,
75  const Tag* file_id_tag )
76 {
77  // At this time there is no support for reading a subset of the file
78  if( subset_list )
79  {
80  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for NASTRAN" );
81  }
82 
83  nodeIdMap.clear();
84  elemIdMap.clear();
85 
86  bool debug = false;
87  if( debug ) std::cout << "begin ReadNASTRAN::load_file" << std::endl;
88  ErrorCode result;
89 
90  // Count the entities of each type in the file. This is used to allocate the node array.
91  int entity_count[MBMAXTYPE];
92  for( int i = 0; i < MBMAXTYPE; i++ )
93  entity_count[i] = 0;
94 
95  /* Determine the line_format of the first line. Assume that the entire file
96  has the same format. */
97  std::string line;
98  std::ifstream file( filename );
99  if( !getline( file, line ) ) return MB_FILE_DOES_NOT_EXIST;
100  line_format format;
101  result = determine_line_format( line, format );
102  if( MB_SUCCESS != result ) return result;
103 
104  /* Count the number of each entity in the file. This allows one to allocate
105  a sequential array of vertex handles. */
106  while( !file.eof() )
107  {
108  // Cut the line into fields as determined by the line format.
109  // Use a vector to allow for an unknown number of tokens (continue lines).
110  // Continue lines are not implemented.
111  std::vector< std::string > tokens;
112  tokens.reserve( 10 ); // assume 10 fields to avoid extra vector resizing
113  result = tokenize_line( line, format, tokens );
114  if( MB_SUCCESS != result ) return result;
115 
116  // Process the tokens of the line. The first token describes the entity type.
117  EntityType type;
118  result = determine_entity_type( ( tokens.empty() ) ? "" : tokens.front(), type );
119  if( MB_SUCCESS != result ) return result;
120  entity_count[type]++;
121  getline( file, line );
122  }
123 
124  if( debug )
125  {
126  for( int i = 0; i < MBMAXTYPE; i++ )
127  {
128  std::cout << "entity_count[" << i << "]=" << entity_count[i] << std::endl;
129  }
130  }
131 
132  // Keep list of material sets
133  std::vector< Range > materials;
134 
135  // Now that the number of vertices is known, create the vertices.
136  EntityHandle start_vert = 0;
137  std::vector< double* > coord_arrays( 3 );
138  result = readMeshIface->get_node_coords( 3, entity_count[0], MB_START_ID, start_vert, coord_arrays );
139  if( MB_SUCCESS != result ) return result;
140  if( 0 == start_vert ) return MB_FAILURE; // check for NULL
141  int id, vert_index = 0;
142  if( debug ) std::cout << "allocated coord arrays" << std::endl;
143 
144  // Read the file again to create the entities.
145  file.clear(); // Clear eof state from object
146  file.seekg( 0 ); // Rewind file
147  while( !file.eof() )
148  {
149  getline( file, line );
150 
151  // Cut the line into fields as determined by the line format.
152  // Use a vector to allow for an unknown number of tokens (continue lines).
153  // Continue lines are not implemented.
154  std::vector< std::string > tokens;
155  tokens.reserve( 10 ); // assume 10 fields to avoid extra vector resizing
156  result = tokenize_line( line, format, tokens );
157  if( MB_SUCCESS != result ) return result;
158 
159  // Process the tokens of the line. The first token describes the entity type.
160  EntityType type;
161  result = determine_entity_type( tokens.front(), type );
162  if( MB_SUCCESS != result ) return result;
163 
164  // Create the entity.
165  if( MBVERTEX == type )
166  {
167  double* coords[3] = { coord_arrays[0] + vert_index, coord_arrays[1] + vert_index,
168  coord_arrays[2] + vert_index };
169  result = read_node( tokens, debug, coords, id );
170  if( MB_SUCCESS != result ) return result;
171  if( !nodeIdMap.insert( id, start_vert + vert_index, 1 ).second ) return MB_FAILURE; // Duplicate IDs!
172  ++vert_index;
173  }
174  else
175  {
176  result = read_element( tokens, materials, type, debug );
177  if( MB_SUCCESS != result ) return result;
178  }
179  }
180 
181  result = create_materials( materials );
182  if( MB_SUCCESS != result ) return result;
183 
184  result = assign_ids( file_id_tag );
185  if( MB_SUCCESS != result ) return result;
186 
187  file.close();
188  nodeIdMap.clear();
189  elemIdMap.clear();
190  return MB_SUCCESS;
191 }
192 
193 /* Determine the type of NASTRAN line: small field, large field, or free field.
194  small field: each line has 10 fields of 8 characters
195  large field: 1x8, 4x16, 1x8. Field 1 must have an asterisk following the character string
196  free field: each line entry must be separated by a comma
197  Implementation tries to avoid more searches than necessary. */
198 ErrorCode ReadNASTRAN::determine_line_format( const std::string& line, line_format& format )
199 {
200  std::string::size_type found_asterisk = line.find( "*" );
201  if( std::string::npos != found_asterisk )
202  {
203  format = LARGE_FIELD;
204  return MB_SUCCESS;
205  }
206  else
207  {
208  std::string::size_type found_comma = line.find( "," );
209  if( std::string::npos != found_comma )
210  {
211  format = FREE_FIELD;
212  return MB_SUCCESS;
213  }
214  else
215  {
216  format = SMALL_FIELD;
217  return MB_SUCCESS;
218  }
219  }
220 }
221 
222 /* Tokenize the line. Continue-lines have not been implemented. */
223 ErrorCode ReadNASTRAN::tokenize_line( const std::string& line,
224  const line_format format,
225  std::vector< std::string >& tokens )
226 {
227  size_t line_size = line.size();
228  switch( format )
229  {
230  case SMALL_FIELD: {
231  // Expect 10 fields of 8 characters.
232  // The sample file does not have all 10 fields in each line
233  const int field_length = 8;
234  unsigned int n_tokens = line_size / field_length;
235  for( unsigned int i = 0; i < n_tokens; i++ )
236  {
237  tokens.push_back( line.substr( i * field_length, field_length ) );
238  }
239  break;
240  }
241  case LARGE_FIELD:
242  return MB_NOT_IMPLEMENTED;
243  case FREE_FIELD:
244  return MB_NOT_IMPLEMENTED;
245  default:
246  return MB_FAILURE;
247  }
248 
249  return MB_SUCCESS;
250 }
251 
252 ErrorCode ReadNASTRAN::determine_entity_type( const std::string& first_token, EntityType& type )
253 {
254  if( 0 == first_token.compare( "GRID " ) )
255  type = MBVERTEX;
256  else if( 0 == first_token.compare( "CTETRA " ) )
257  type = MBTET;
258  else if( 0 == first_token.compare( "CPENTA " ) )
259  type = MBPRISM;
260  else if( 0 == first_token.compare( "CHEXA " ) )
261  type = MBHEX;
262  else
263  return MB_NOT_IMPLEMENTED;
264 
265  return MB_SUCCESS;
266 }
267 
268 /* Some help from Jason:
269  Nastran floats must contain a decimal point, may contain
270  a leading '-' and may contain an exponent. The 'E' is optional
271  when specifying an exponent. A '-' or '+' at any location other
272  than the first position indicates an exponent. For a positive
273  exponent, either a '+' or an 'E' must be specified. For a
274  negative exponent, the 'E' is option and the '-' is always specified.
275  Examples for the real value 7.0 from mcs2006 quick reference guide:
276  7.0 .7E1 0.7+1 .70+1 7.E+0 70.-1
277 
278  From the test file created in SC/Tetra:
279  GRID 1 03.9804546.9052-15.6008-1
280  has the coordinates: ( 3.980454, 6.9052e-1, 5.6008e-1 )
281  GRID 200005 04.004752-3.985-15.4955-1
282  has the coordinates: ( 4.004752, -3.985e-1, 5.4955e-1 ) */
283 ErrorCode ReadNASTRAN::get_real( const std::string& token, double& real )
284 {
285  std::string significand = token;
286  std::string exponent = "0";
287 
288  // Cut off the first digit because a "-" could be here indicating a negative
289  // number. Instead we are looking for a negative exponent.
290  std::string back_token = token.substr( 1 );
291 
292  // A minus that is not the first digit is always a negative exponent
293  std::string::size_type found_minus = back_token.find( "-" );
294  if( std::string::npos != found_minus )
295  {
296  // separate the significand from the exponent at the "-"
297  exponent = token.substr( found_minus + 1 );
298  significand = token.substr( 0, found_minus + 1 );
299 
300  // If the significand has an "E", remove it
301  if( std::string::npos != significand.find( "E" ) )
302  // Assume the "E" is at the end of the significand.
303  significand = significand.substr( 1, significand.size() - 2 );
304 
305  // If a minus does not exist past the 1st digit, but an "E" or "+" does, then
306  // it is a positive exponent. First look for an "E",
307  }
308  else
309  {
310  std::string::size_type found_E = token.find( "E" );
311  if( std::string::npos != found_E )
312  {
313  significand = token.substr( 0, found_E - 1 );
314  exponent = token.substr( found_E + 1 );
315  // If there is a "+" on the exponent, cut it off
316  std::size_t found_plus = exponent.find( "+" );
317  if( std::string::npos != found_plus )
318  {
319  exponent = exponent.substr( found_plus + 1 );
320  }
321  }
322  else
323  {
324  // If there is a "+" on the exponent, cut it off
325  std::size_t found_plus = token.find( "+" );
326  if( std::string::npos != found_plus )
327  {
328  significand = token.substr( 0, found_plus - 1 );
329  exponent = token.substr( found_plus + 1 );
330  }
331  }
332  }
333 
334  // Now assemble the real number
335  double signi = atof( significand.c_str() );
336  double expon = atof( exponent.c_str() );
337 
338  if( HUGE_VAL == signi || HUGE_VAL == expon ) return MB_FAILURE;
339 
340  real = signi * pow( 10, expon );
341 
342  return MB_SUCCESS;
343 }
344 
345 /* It has been determined that this line is a vertex. Read the rest of
346  the line and create the vertex. */
347 ErrorCode ReadNASTRAN::read_node( const std::vector< std::string >& tokens,
348  const bool debug,
349  double* coords[3],
350  int& id )
351 {
352  // Read the node's id (unique)
353  ErrorCode result;
354  id = atoi( tokens[1].c_str() );
355 
356  // Read the node's coordinate system number
357  // "0" or blank refers to the basic coordinate system.
358  int coord_system = atoi( tokens[2].c_str() );
359  if( 0 != coord_system )
360  {
361  std::cerr << "ReadNASTRAN: alternative coordinate systems not implemented" << std::endl;
362  return MB_NOT_IMPLEMENTED;
363  }
364 
365  // Read the coordinates
366  for( unsigned int i = 0; i < 3; i++ )
367  {
368  result = get_real( tokens[i + 3], *coords[i] );
369  if( MB_SUCCESS != result ) return result;
370  if( debug ) std::cout << "read_node: coords[" << i << "]=" << coords[i] << std::endl;
371  }
372 
373  return MB_SUCCESS;
374 }
375 
376 /* The type of element has already been identified. Read the rest of the
377  line and create the element. Assume that all of the nodes have already
378  been created. */
379 ErrorCode ReadNASTRAN::read_element( const std::vector< std::string >& tokens,
380  std::vector< Range >& materials,
381  const EntityType element_type,
382  const bool /*debug*/ )
383 {
384  // Read the element's id (unique) and material set
385  ErrorCode result;
386  int id = atoi( tokens[1].c_str() );
387  int material = atoi( tokens[2].c_str() );
388 
389  // Resize materials list if necessary. This code is somewhat complicated
390  // so as to avoid copying of Ranges
391  if( material >= (int)materials.size() )
392  {
393  if( (int)materials.capacity() < material )
394  materials.resize( material + 1 );
395  else
396  {
397  std::vector< Range > new_mat( material + 1 );
398  for( size_t i = 0; i < materials.size(); ++i )
399  new_mat[i].swap( materials[i] );
400  materials.swap( new_mat );
401  }
402  }
403 
404  // The size of the connectivity array depends on the element type
405  int n_conn = CN::VerticesPerEntity( element_type );
406  EntityHandle conn_verts[27];
407  assert( n_conn <= (int)( sizeof( conn_verts ) / sizeof( EntityHandle ) ) );
408 
409  // Read the connected node ids from the file
410  for( int i = 0; i < n_conn; i++ )
411  {
412  int n = atoi( tokens[3 + i].c_str() );
413  conn_verts[i] = nodeIdMap.find( n );
414  if( !conn_verts[i] ) // invalid vertex id
415  return MB_FAILURE;
416  }
417 
418  // Create the element and set the global id from the NASTRAN file
419  EntityHandle element;
420  result = MBI->create_element( element_type, conn_verts, n_conn, element );
421  if( MB_SUCCESS != result ) return result;
422  elemIdMap.insert( id, element, 1 );
423 
424  materials[material].insert( element );
425  return MB_SUCCESS;
426 }
427 
428 ErrorCode ReadNASTRAN::create_materials( const std::vector< Range >& materials )
429 {
430  ErrorCode result;
431  Tag material_tag;
432  int negone = -1;
434  &negone );
435  if( MB_SUCCESS != result ) return result;
436 
437  for( size_t i = 0; i < materials.size(); ++i )
438  {
439  if( materials[i].empty() ) continue;
440 
441  // Merge with existing or create new? Original code effectively
442  // created new by only merging with existing in current file set,
443  // so do the same here. - j.kraftcheck
444 
445  EntityHandle handle;
446  result = MBI->create_meshset( MESHSET_SET, handle );
447  if( MB_SUCCESS != result ) return result;
448 
449  result = MBI->add_entities( handle, materials[i] );
450  if( MB_SUCCESS != result ) return result;
451 
452  int id = i;
453  result = MBI->tag_set_data( material_tag, &handle, 1, &id );
454  if( MB_SUCCESS != result ) return result;
455  }
456 
457  return MB_SUCCESS;
458 }
459 
460 ErrorCode ReadNASTRAN::assign_ids( const Tag* file_id_tag )
461 {
462  // Create tag
463  ErrorCode result;
464  Tag id_tag = MBI->globalId_tag();
465 
467  for( int t = 0; t < 2; ++t )
468  {
470  for( i = fileIdMap.begin(); i != fileIdMap.end(); ++i )
471  {
472  Range range( i->value, i->value + i->count - 1 );
473 
474  result = readMeshIface->assign_ids( id_tag, range, i->begin );
475  if( MB_SUCCESS != result ) return result;
476 
477  if( file_id_tag && *file_id_tag != id_tag )
478  {
479  result = readMeshIface->assign_ids( *file_id_tag, range, i->begin );
480  if( MB_SUCCESS != result ) return result;
481  }
482  }
483  }
484 
485  return MB_SUCCESS;
486 }
487 
488 } // namespace moab