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
/**
 * 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.
 *
 */

#ifdef WIN32
#ifdef _DEBUG
// turn off warnings that say they debugging identifier has been truncated
// this warning comes up when using some STL containers
#pragma warning( disable : 4786 )
#endif
#endif

#include "WriteAns.hpp"

#include <utility>
#include <algorithm>
#include <ctime>
#include <string>
#include <vector>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <iomanip>

#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include <cassert>
#include "Internals.hpp"
#include "ExoIIUtil.hpp"
#include "MBTagConventions.hpp"

#define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )

namespace moab
{

WriterIface* WriteAns::factory( Interface* iface )
{
    return new WriteAns( iface );
}

WriteAns::WriteAns( Interface* impl ) : mbImpl( impl ), mCurrentMeshHandle( 0 ), mGlobalIdTag( 0 ), mMatSetIdTag( 0 )
{
    assert( impl != NULL );

    // impl->query_interface( mWriteIface );

    // initialize in case tag_get_handle fails below
    //! get and cache predefined tag handles
    const int negone = -1;
    impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                          &negone );

    impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                          &negone );

    impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                          &negone );
}

WriteAns::~WriteAns()
{
    // mbImpl->release_interface(mWriteIface);
}

ErrorCode WriteAns::write_file( const char* file_name,
                                const bool /* overwrite (commented out to remove warning) */,
                                const FileOptions&,
                                const EntityHandle* ent_handles,
                                const int num_sets,
                                const std::vector< std::string >&,
                                const Tag*,
                                int,
                                int )
{
    assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );

    ErrorCode result;

    // set SOLID45 element type to #60000, hope nobody has already...
    const char* ETSolid45 = "60045";
    const char* ETSolid92 = "60042";
    const char* ETSolid95 = "60095";

    // set Material id # to be used as default for all elements
    // will need to be subsequently reassigned inside ANSYS
    // Can, although have not, declare similar defaults for other attributes
    const char* MATDefault = "1";

    // create file streams for writing
    std::ofstream node_file;
    std::ofstream elem_file;
    std::ofstream ans_file;

    // get base filename from filename.ans
    std::string temp_string;
    std::string base_string;
    base_string.assign( file_name );
    base_string.replace( base_string.find_last_of( ".ans" ) - 3, 4, "" );

    // open node file for writing
    temp_string = base_string + ".node";
    node_file.open( temp_string.c_str() );
    node_file.setf( std::ios::scientific, std::ios::floatfield );
    node_file.precision( 13 );

    // open elem file for writing
    temp_string = base_string + ".elem";
    elem_file.open( temp_string.c_str() );

    // open ans file for writing
    ans_file.open( file_name );
    ans_file << "/prep7" << std::endl;

    // gather single output set
    EntityHandle output_set = 0;
    if( ent_handles && num_sets > 0 )
    {
        for( int i = 0; i < num_sets; i++ )
        {
            // from template, maybe can be removed
            result = mbImpl->unite_meshset( output_set, ent_handles[i] );
            if( result != MB_SUCCESS ) return result;
        }
    }

    // search for all nodes
    Range node_range;
    result = mbImpl->get_entities_by_type( output_set, MBVERTEX, node_range, true );
    if( result != MB_SUCCESS ) return result;

    // Commented out until Seg Fault taken care of in gather_nodes...
    // get any missing nodes which are needed for elements
    // Range all_ent_range,missing_range;
    // result=mbImpl->get_entities_by_handle(output_set,all_ent_range,true);
    // if(result !=MB_SUCCESS) return result;
    // result=mWriteIface->gather_nodes_from_elements(all_ent_range,0,missing_range);
    // node_range.merge(missing_range);

    // write the nodes
    double coord[3];
    for( Range::iterator it = node_range.begin(); it != node_range.end(); ++it )
    {
        EntityHandle node_handle = *it;

        result = mbImpl->get_coords( &node_handle, 1, coord );
        if( result != MB_SUCCESS ) return result;

        node_file.width( 8 );
        node_file << mbImpl->id_from_handle( node_handle );
        node_file.width( 20 );
        node_file << coord[0];
        node_file.width( 20 );
        node_file << coord[1];
        node_file.width( 20 );
        node_file << coord[2] << std::endl;
    }

    // update header to load nodes
    ans_file << "nread," << base_string << ",node" << std::endl;

    // search for all node sets (Dirichlet Sets)
    Range node_mesh_sets;
    int ns_id;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, node_mesh_sets );
    if( result != MB_SUCCESS ) return result;

    for( Range::iterator ns_it = node_mesh_sets.begin(); ns_it != node_mesh_sets.end(); ++ns_it )
    {
        result = mbImpl->tag_get_data( mDirichletSetTag, &( *ns_it ), 1, &ns_id );
        if( result != MB_SUCCESS ) return result;
        std::vector< EntityHandle > node_vector;
        result = mbImpl->get_entities_by_handle( *ns_it, node_vector, true );
        if( result != MB_SUCCESS ) return result;
        // for every nodeset found, cycle through nodes in set:
        for( std::vector< EntityHandle >::iterator node_it = node_vector.begin(); node_it != node_vector.end();
             ++node_it )
        {
            int ns_node_id = mbImpl->id_from_handle( *node_it );
            if( node_it == node_vector.begin() )
            {
                // select first node in new list
                ans_file << "nsel,s,node,," << std::setw( 8 ) << ns_node_id << std::endl;
            }
            else
            {
                // append node to list
                ans_file << "nsel,a,node,," << std::setw( 8 ) << ns_node_id << std::endl;
            }
        }
        // create NS(#) nodeset
        ans_file << "cm,NS" << ns_id << ",node" << std::endl;
    }

    // ANSYS Element format:
    // I, J, K, L, M, N, O, P,etc... MAT, TYPE, REAL, SECNUM, ESYS, IEL
    // I-P are nodes of element
    // MAT = material number
    // TYPE = Element type number
    // REAL = Real constant set number
    // SECNUM = section attribute number
    // ESYS = coordinate system for nodes
    // IEL = element # (unique?)
    // For all nodes past 8, write on second line

    // Write all MBTET elements
    Range tet_range;
    result = mbImpl->get_entities_by_type( output_set, MBTET, tet_range, true );
    if( result != MB_SUCCESS ) return result;
    for( Range::iterator elem_it = tet_range.begin(); elem_it != tet_range.end(); ++elem_it )
    {
        EntityHandle elem_handle = *elem_it;
        int elem_id              = mbImpl->id_from_handle( elem_handle );
        std::vector< EntityHandle > conn;
        result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
        if( result != MB_SUCCESS ) return result;
        // make sure 4 or 10 node tet
        if( conn.size() != 4 && conn.size() != 10 )
        {
            std::cout << "Support not added for element type. \n";
            return MB_FAILURE;
        }
        // write information for 4 node tet
        if( conn.size() == 4 )
        {
            elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
            elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[2];
            elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];
            elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];

            elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
            elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
            elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
            elem_file << std::endl;
        }

        // write information for 10 node tet
        if( conn.size() == 10 )
        {
            elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
            elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
            elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
            elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];

            elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid92;
            elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
            elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
            elem_file << std::endl;

            elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[9];
            elem_file << std::endl;
        }
    }

    // Write all MBHEX elements
    Range hex_range;
    result = mbImpl->get_entities_by_type( output_set, MBHEX, hex_range, true );
    if( result != MB_SUCCESS ) return result;
    for( Range::iterator elem_it = hex_range.begin(); elem_it != hex_range.end(); ++elem_it )
    {
        EntityHandle elem_handle = *elem_it;
        int elem_id              = mbImpl->id_from_handle( elem_handle );
        std::vector< EntityHandle > conn;
        result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
        if( result != MB_SUCCESS ) return result;
        // make sure supported hex type
        if( conn.size() != 8 && conn.size() != 20 )
        {
            std::cout << "Support not added for element type. \n";
            return MB_FAILURE;
        }

        // write information for 8 node hex
        if( conn.size() == 8 )
        {
            elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
            elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
            elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
            elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];

            elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
            elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
            elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
            elem_file << std::endl;
        }

        // write information for 20 node hex
        if( conn.size() == 20 )
        {

            elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
            elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[0];
            elem_file << std::setw( 8 ) << conn[7] << std::setw( 8 ) << conn[6];
            elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];

            elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid95;
            elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
            elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
            elem_file << std::endl;

            elem_file << std::setw( 8 ) << conn[16] << std::setw( 8 ) << conn[13];
            elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[12];
            elem_file << std::setw( 8 ) << conn[18] << std::setw( 8 ) << conn[14];
            elem_file << std::setw( 8 ) << conn[10] << std::setw( 8 ) << conn[15];
            elem_file << std::setw( 8 ) << conn[19] << std::setw( 8 ) << conn[17];
            elem_file << std::setw( 8 ) << conn[9] << std::setw( 8 ) << conn[11];
            elem_file << std::endl;
        }
    }
    // Write all MBPRISM elements
    Range prism_range;
    result = mbImpl->get_entities_by_type( output_set, MBPRISM, prism_range, true );
    if( result != MB_SUCCESS ) return result;
    for( Range::iterator elem_it = prism_range.begin(); elem_it != prism_range.end(); ++elem_it )
    {
        EntityHandle elem_handle = *elem_it;
        int elem_id              = mbImpl->id_from_handle( elem_handle );
        std::vector< EntityHandle > conn;
        result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
        if( result != MB_SUCCESS ) return result;
        // make sure supported prism type
        if( conn.size() != 6 )<--- Assuming that condition 'conn.size()!=6' is not redundant
        {
            std::cout << "Support not added for element type. \n";
            return MB_FAILURE;
        }

        // write information for 6 node prism
        if( conn.size() == 6 )<--- Condition 'conn.size()==6' is always true
        {
            elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[3];
            elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[4];
            elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[2];
            elem_file << std::setw( 8 ) << conn[5] << std::setw( 8 ) << conn[5];

            elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
            elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
            elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
            elem_file << std::endl;
        }
    }

    // create element types (for now writes all, even if not used)
    ans_file << "et," << ETSolid45 << ",SOLID45" << std::endl;
    ans_file << "et," << ETSolid92 << ",SOLID92" << std::endl;
    ans_file << "et," << ETSolid95 << ",SOLID95" << std::endl;

    // xxx pyramids, other elements later...

    // write header to load elements
    ans_file << "eread," << base_string << ",elem" << std::endl;

    // search for all side sets (Neumann)
    Range side_mesh_sets;
    int ss_id;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, side_mesh_sets );
    if( result != MB_SUCCESS ) return result;
    // cycle through all sets found
    for( Range::iterator ss_it = side_mesh_sets.begin(); ss_it != side_mesh_sets.end(); ++ss_it )
    {
        result = mbImpl->tag_get_data( mNeumannSetTag, &( *ss_it ), 1, &ss_id );
        if( result != MB_SUCCESS ) return result;
        std::vector< EntityHandle > elem_vector;
        result = mbImpl->get_entities_by_handle( *ss_it, elem_vector, true );
        if( result != MB_SUCCESS ) return result;

        // cycle through elements in current side set
        for( std::vector< EntityHandle >::iterator elem_it = elem_vector.begin(); elem_it != elem_vector.end();
             ++elem_it )
        {
            EntityHandle elem_handle = *elem_it;

            // instead of selecting current element in set, select its nodes...
            std::vector< EntityHandle > conn;
            result = mbImpl->get_connectivity( &elem_handle, 1, conn );
            if( result != MB_SUCCESS ) return result;
            if( elem_it == elem_vector.begin() )
            {
                ans_file << "nsel,s,node,," << std::setw( 8 ) << conn[0] << std::endl;
                for( unsigned int i = 1; i < conn.size(); i++ )
                {
                    ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
                }
            }
            else
            {
                for( unsigned int i = 0; i < conn.size(); i++ )
                {
                    ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
                }
            }
        }
        // create SS(#) node set
        ans_file << "cm,SS" << ss_id << ",node" << std::endl;
    }

    // Gather all element blocks
    Range matset;
    int mat_id;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, matset );
    if( result != MB_SUCCESS ) return result;
    // cycle through all elem blocks
    for( Range::iterator mat_it = matset.begin(); mat_it != matset.end(); ++mat_it )
    {
        EntityHandle matset_handle = *mat_it;
        result                     = mbImpl->tag_get_data( mMaterialSetTag, &matset_handle, 1, &mat_id );
        if( result != MB_SUCCESS ) return result;
        std::vector< EntityHandle > mat_vector;
        result = mbImpl->get_entities_by_handle( *mat_it, mat_vector, true );
        if( result != MB_SUCCESS ) return result;
        // cycle through elements in current mat set
        for( std::vector< EntityHandle >::iterator elem_it = mat_vector.begin(); elem_it != mat_vector.end();
             ++elem_it )
        {
            EntityHandle elem_handle = *elem_it;
            int elem_id              = mbImpl->id_from_handle( elem_handle );
            if( elem_it == mat_vector.begin() )
            {
                ans_file << "esel,s,elem,," << std::setw( 8 ) << elem_id << std::endl;
            }
            else
            {
                ans_file << "esel,a,elem,," << std::setw( 8 ) << elem_id << std::endl;
            }
        }
        // for each matset, write block command
        ans_file << "cm,EB" << mat_id << ",elem" << std::endl;
    }

    // close all file streams
    node_file.close();
    elem_file.close();
    ans_file.close();

    return MB_SUCCESS;
}

}  // namespace moab