1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

#include "ReadNASTRAN.hpp"

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>
#include <cassert>
#include <cmath>

#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "Internals.hpp"  // For MB_START_ID
#include "moab/Range.hpp"
#include "moab/FileOptions.hpp"
#include "FileTokenizer.hpp"
#include "MBTagConventions.hpp"
#include "moab/CN.hpp"

namespace moab
{

ReaderIface* ReadNASTRAN::factory( Interface* iface )
{
    return new ReadNASTRAN( iface );
}

// Constructor
ReadNASTRAN::ReadNASTRAN( Interface* impl ) : MBI( impl )
{
    assert( NULL != impl );
    MBI->query_interface( readMeshIface );
    assert( NULL != readMeshIface );
}

// Destructor
ReadNASTRAN::~ReadNASTRAN()
{
    if( readMeshIface )
    {
        MBI->release_interface( readMeshIface );
        readMeshIface = 0;
    }
}

ErrorCode ReadNASTRAN::read_tag_values( const char* /*file_name*/,
                                        const char* /*tag_name*/,
                                        const FileOptions& /*opts*/,
                                        std::vector< int >& /*tag_values_out*/,
                                        const SubsetList* /*subset_list*/ )
{
    return MB_NOT_IMPLEMENTED;
}

// Load the file as called by the Interface function
ErrorCode ReadNASTRAN::load_file( const char* filename,
                                  const EntityHandle* /* file_set */,
                                  const FileOptions& /* opts */,
                                  const ReaderIface::SubsetList* subset_list,
                                  const Tag* file_id_tag )
{
    // At this time there is no support for reading a subset of the file
    if( subset_list )
    {
        MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for NASTRAN" );
    }

    nodeIdMap.clear();
    elemIdMap.clear();

    bool debug = false;<--- Assignment 'debug=false', assigned value is 0
    if( debug ) std::cout << "begin ReadNASTRAN::load_file" << std::endl;<--- Condition 'debug' is always false
    ErrorCode result;

    // Count the entities of each type in the file. This is used to allocate the node array.
    int entity_count[MBMAXTYPE];
    for( int i = 0; i < MBMAXTYPE; i++ )
        entity_count[i] = 0;

    /* Determine the line_format of the first line. Assume that the entire file
       has the same format. */
    std::string line;
    std::ifstream file( filename );
    if( !getline( file, line ) ) return MB_FILE_DOES_NOT_EXIST;
    line_format format;
    result = determine_line_format( line, format );
    if( MB_SUCCESS != result ) return result;

    /* Count the number of each entity in the file. This allows one to allocate
       a sequential array of vertex handles. */
    while( !file.eof() )
    {
        // Cut the line into fields as determined by the line format.
        // Use a vector to allow for an unknown number of tokens (continue lines).
        // Continue lines are not implemented.
        std::vector< std::string > tokens;
        tokens.reserve( 10 );  // assume 10 fields to avoid extra vector resizing
        result = tokenize_line( line, format, tokens );
        if( MB_SUCCESS != result ) return result;

        // Process the tokens of the line. The first token describes the entity type.
        EntityType type;
        result = determine_entity_type( ( tokens.empty() ) ? "" : tokens.front(), type );
        if( MB_SUCCESS != result ) return result;
        entity_count[type]++;
        getline( file, line );
    }

    if( debug )
    {
        for( int i = 0; i < MBMAXTYPE; i++ )
        {
            std::cout << "entity_count[" << i << "]=" << entity_count[i] << std::endl;
        }
    }

    // Keep list of material sets
    std::vector< Range > materials;

    // Now that the number of vertices is known, create the vertices.
    EntityHandle start_vert = 0;
    std::vector< double* > coord_arrays( 3 );
    result = readMeshIface->get_node_coords( 3, entity_count[0], MB_START_ID, start_vert, coord_arrays );
    if( MB_SUCCESS != result ) return result;
    if( 0 == start_vert ) return MB_FAILURE;  // check for NULL
    int id, vert_index = 0;
    if( debug ) std::cout << "allocated coord arrays" << std::endl;

    // Read the file again to create the entities.
    file.clear();     // Clear eof state from object
    file.seekg( 0 );  // Rewind file
    while( !file.eof() )
    {
        getline( file, line );

        // Cut the line into fields as determined by the line format.
        // Use a vector to allow for an unknown number of tokens (continue lines).
        // Continue lines are not implemented.
        std::vector< std::string > tokens;
        tokens.reserve( 10 );  // assume 10 fields to avoid extra vector resizing
        result = tokenize_line( line, format, tokens );
        if( MB_SUCCESS != result ) return result;

        // Process the tokens of the line. The first token describes the entity type.
        EntityType type;
        result = determine_entity_type( tokens.front(), type );
        if( MB_SUCCESS != result ) return result;

        // Create the entity.
        if( MBVERTEX == type )
        {
            double* coords[3] = { coord_arrays[0] + vert_index, coord_arrays[1] + vert_index,
                                  coord_arrays[2] + vert_index };
            result            = read_node( tokens, debug, coords, id );
            if( MB_SUCCESS != result ) return result;
            if( !nodeIdMap.insert( id, start_vert + vert_index, 1 ).second ) return MB_FAILURE;  // Duplicate IDs!
            ++vert_index;
        }
        else
        {
            result = read_element( tokens, materials, type, debug );
            if( MB_SUCCESS != result ) return result;
        }
    }

    result = create_materials( materials );
    if( MB_SUCCESS != result ) return result;

    result = assign_ids( file_id_tag );
    if( MB_SUCCESS != result ) return result;

    file.close();
    nodeIdMap.clear();
    elemIdMap.clear();
    return MB_SUCCESS;
}

/* Determine the type of NASTRAN line: small field, large field, or free field.
   small field: each line has 10 fields of 8 characters
   large field: 1x8, 4x16, 1x8. Field 1 must have an asterisk following the character string
   free field: each line entry must be separated by a comma
   Implementation tries to avoid more searches than necessary. */
ErrorCode ReadNASTRAN::determine_line_format( const std::string& line, line_format& format )
{
    std::string::size_type found_asterisk = line.find( "*" );
    if( std::string::npos != found_asterisk )
    {
        format = LARGE_FIELD;
        return MB_SUCCESS;
    }
    else
    {
        std::string::size_type found_comma = line.find( "," );
        if( std::string::npos != found_comma )
        {
            format = FREE_FIELD;
            return MB_SUCCESS;
        }
        else
        {
            format = SMALL_FIELD;
            return MB_SUCCESS;
        }
    }
}

/* Tokenize the line. Continue-lines have not been implemented. */
ErrorCode ReadNASTRAN::tokenize_line( const std::string& line,
                                      const line_format format,
                                      std::vector< std::string >& tokens )
{
    size_t line_size = line.size();
    switch( format )
    {
        case SMALL_FIELD: {
            // Expect 10 fields of 8 characters.
            // The sample file does not have all 10 fields in each line
            const int field_length = 8;
            unsigned int n_tokens  = line_size / field_length;
            for( unsigned int i = 0; i < n_tokens; i++ )
            {
                tokens.push_back( line.substr( i * field_length, field_length ) );
            }
            break;
        }
        case LARGE_FIELD:
            return MB_NOT_IMPLEMENTED;
        case FREE_FIELD:
            return MB_NOT_IMPLEMENTED;
        default:
            return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode ReadNASTRAN::determine_entity_type( const std::string& first_token, EntityType& type )
{
    if( 0 == first_token.compare( "GRID    " ) )
        type = MBVERTEX;
    else if( 0 == first_token.compare( "CTETRA  " ) )
        type = MBTET;
    else if( 0 == first_token.compare( "CPENTA  " ) )
        type = MBPRISM;
    else if( 0 == first_token.compare( "CHEXA   " ) )
        type = MBHEX;
    else
        return MB_NOT_IMPLEMENTED;

    return MB_SUCCESS;
}

/* Some help from Jason:
   Nastran floats must contain a decimal point, may contain
   a leading '-' and may contain an exponent. The 'E' is optional
   when specifying an exponent.  A '-' or '+' at any location other
   than the first position indicates an exponent.  For a positive
   exponent, either a '+' or an 'E' must be specified.  For a
   negative exponent, the 'E' is option and the '-' is always specified.
   Examples for the real value 7.0 from mcs2006 quick reference guide:
   7.0  .7E1  0.7+1  .70+1  7.E+0  70.-1

   From the test file created in SC/Tetra:
   GRID           1       03.9804546.9052-15.6008-1
   has the coordinates: ( 3.980454, 6.9052e-1, 5.6008e-1 )
   GRID      200005       04.004752-3.985-15.4955-1
   has the coordinates: ( 4.004752, -3.985e-1, 5.4955e-1 ) */
ErrorCode ReadNASTRAN::get_real( const std::string& token, double& real )
{
    std::string significand = token;
    std::string exponent    = "0";

    // Cut off the first digit because a "-" could be here indicating a negative
    // number. Instead we are looking for a negative exponent.
    std::string back_token = token.substr( 1 );

    // A minus that is not the first digit is always a negative exponent
    std::string::size_type found_minus = back_token.find( "-" );
    if( std::string::npos != found_minus )
    {
        // separate the significand from the exponent at the "-"
        exponent    = token.substr( found_minus + 1 );
        significand = token.substr( 0, found_minus + 1 );

        // If the significand has an "E", remove it
        if( std::string::npos != significand.find( "E" ) )
            // Assume the "E" is at the end of the significand.
            significand = significand.substr( 1, significand.size() - 2 );

        // If a minus does not exist past the 1st digit, but an "E" or "+" does, then
        // it is a positive exponent. First look for an "E",
    }
    else
    {
        std::string::size_type found_E = token.find( "E" );
        if( std::string::npos != found_E )
        {
            significand = token.substr( 0, found_E - 1 );
            exponent    = token.substr( found_E + 1 );
            // If there is a "+" on the exponent, cut it off
            std::size_t found_plus = exponent.find( "+" );
            if( std::string::npos != found_plus )
            {
                exponent = exponent.substr( found_plus + 1 );
            }
        }
        else
        {
            // If there is a "+" on the exponent, cut it off
            std::size_t found_plus = token.find( "+" );
            if( std::string::npos != found_plus )
            {
                significand = token.substr( 0, found_plus - 1 );
                exponent    = token.substr( found_plus + 1 );
            }
        }
    }

    // Now assemble the real number
    double signi = atof( significand.c_str() );
    double expon = atof( exponent.c_str() );

    if( HUGE_VAL == signi || HUGE_VAL == expon ) return MB_FAILURE;

    real = signi * pow( 10, expon );

    return MB_SUCCESS;
}

/* It has been determined that this line is a vertex. Read the rest of
   the line and create the vertex. */
ErrorCode ReadNASTRAN::read_node( const std::vector< std::string >& tokens,
                                  const bool debug,
                                  double* coords[3],
                                  int& id )
{
    // Read the node's id (unique)
    ErrorCode result;
    id = atoi( tokens[1].c_str() );

    // Read the node's coordinate system number
    // "0" or blank refers to the basic coordinate system.
    int coord_system = atoi( tokens[2].c_str() );
    if( 0 != coord_system )
    {
        std::cerr << "ReadNASTRAN: alternative coordinate systems not implemented" << std::endl;
        return MB_NOT_IMPLEMENTED;
    }

    // Read the coordinates
    for( unsigned int i = 0; i < 3; i++ )
    {
        result = get_real( tokens[i + 3], *coords[i] );
        if( MB_SUCCESS != result ) return result;
        if( debug ) std::cout << "read_node: coords[" << i << "]=" << coords[i] << std::endl;
    }

    return MB_SUCCESS;
}

/* The type of element has already been identified. Read the rest of the
   line and create the element. Assume that all of the nodes have already
   been created. */
ErrorCode ReadNASTRAN::read_element( const std::vector< std::string >& tokens,
                                     std::vector< Range >& materials,
                                     const EntityType element_type,
                                     const bool /*debug*/ )
{
    // Read the element's id (unique) and material set
    ErrorCode result;
    int id       = atoi( tokens[1].c_str() );
    int material = atoi( tokens[2].c_str() );

    // Resize materials list if necessary. This code is somewhat complicated
    // so as to avoid copying of Ranges
    if( material >= (int)materials.size() )
    {
        if( (int)materials.capacity() < material )
            materials.resize( material + 1 );
        else
        {
            std::vector< Range > new_mat( material + 1 );
            for( size_t i = 0; i < materials.size(); ++i )
                new_mat[i].swap( materials[i] );
            materials.swap( new_mat );
        }
    }

    // The size of the connectivity array depends on the element type
    int n_conn = CN::VerticesPerEntity( element_type );
    EntityHandle conn_verts[27];
    assert( n_conn <= (int)( sizeof( conn_verts ) / sizeof( EntityHandle ) ) );

    // Read the connected node ids from the file
    for( int i = 0; i < n_conn; i++ )
    {
        int n         = atoi( tokens[3 + i].c_str() );
        conn_verts[i] = nodeIdMap.find( n );
        if( !conn_verts[i] )  // invalid vertex id
            return MB_FAILURE;
    }

    // Create the element and set the global id from the NASTRAN file
    EntityHandle element;
    result = MBI->create_element( element_type, conn_verts, n_conn, element );
    if( MB_SUCCESS != result ) return result;
    elemIdMap.insert( id, element, 1 );

    materials[material].insert( element );
    return MB_SUCCESS;
}

ErrorCode ReadNASTRAN::create_materials( const std::vector< Range >& materials )
{
    ErrorCode result;
    Tag material_tag;
    int negone = -1;
    result = MBI->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, material_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
                                  &negone );
    if( MB_SUCCESS != result ) return result;

    for( size_t i = 0; i < materials.size(); ++i )
    {
        if( materials[i].empty() ) continue;

        // Merge with existing or create new?  Original code effectively
        // created new by only merging with existing in current file set,
        // so do the same here. - j.kraftcheck

        EntityHandle handle;
        result = MBI->create_meshset( MESHSET_SET, handle );
        if( MB_SUCCESS != result ) return result;

        result = MBI->add_entities( handle, materials[i] );
        if( MB_SUCCESS != result ) return result;

        int id = i;
        result = MBI->tag_set_data( material_tag, &handle, 1, &id );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

ErrorCode ReadNASTRAN::assign_ids( const Tag* file_id_tag )
{
    // Create tag
    ErrorCode result;
    Tag id_tag = MBI->globalId_tag();

    RangeMap< int, EntityHandle >::iterator i;
    for( int t = 0; t < 2; ++t )
    {
        RangeMap< int, EntityHandle >& fileIdMap = t ? elemIdMap : nodeIdMap;
        for( i = fileIdMap.begin(); i != fileIdMap.end(); ++i )
        {
            Range range( i->value, i->value + i->count - 1 );

            result = readMeshIface->assign_ids( id_tag, range, i->begin );
            if( MB_SUCCESS != result ) return result;

            if( file_id_tag && *file_id_tag != id_tag )
            {
                result = readMeshIface->assign_ids( *file_id_tag, range, i->begin );
                if( MB_SUCCESS != result ) return result;
            }
        }
    }

    return MB_SUCCESS;
}

}  // namespace moab