MOAB: Mesh Oriented datABase  (version 5.2.1)
WriteAns.cpp
Go to the documentation of this file.
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 #ifdef WIN32
00017 #ifdef _DEBUG
00018 // turn off warnings that say they debugging identifier has been truncated
00019 // this warning comes up when using some STL containers
00020 #pragma warning( disable : 4786 )
00021 #endif
00022 #endif
00023 
00024 #include "WriteAns.hpp"
00025 
00026 #include <utility>
00027 #include <algorithm>
00028 #include <time.h>
00029 #include <string>
00030 #include <vector>
00031 #include <stdio.h>
00032 #include <iostream>
00033 #include <fstream>
00034 #include <iomanip>
00035 
00036 #include "moab/Interface.hpp"
00037 #include "moab/Range.hpp"
00038 #include "assert.h"
00039 #include "Internals.hpp"
00040 #include "ExoIIUtil.hpp"
00041 #include "MBTagConventions.hpp"
00042 
00043 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00044 
00045 namespace moab
00046 {
00047 
00048 WriterIface* WriteAns::factory( Interface* iface )
00049 {
00050     return new WriteAns( iface );
00051 }
00052 
00053 WriteAns::WriteAns( Interface* impl ) : mbImpl( impl ), mCurrentMeshHandle( 0 ), mGlobalIdTag( 0 ), mMatSetIdTag( 0 )
00054 {
00055     assert( impl != NULL );
00056 
00057     // impl->query_interface( mWriteIface );
00058 
00059     // initialize in case tag_get_handle fails below
00060     //! get and cache predefined tag handles
00061     const int negone = -1;
00062     impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00063                           &negone );
00064 
00065     impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00066                           &negone );
00067 
00068     impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00069                           &negone );
00070 }
00071 
00072 WriteAns::~WriteAns()
00073 {
00074     // mbImpl->release_interface(mWriteIface);
00075 }
00076 
00077 ErrorCode WriteAns::write_file( const char* file_name, const bool /* overwrite (commented out to remove warning) */,
00078                                 const FileOptions&, const EntityHandle* ent_handles, const int num_sets,
00079                                 const std::vector< std::string >&, const Tag*, int, int )
00080 {
00081     assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
00082 
00083     ErrorCode result;
00084 
00085     // set SOLID45 element type to #60000, hope nobody has already...
00086     const char* ETSolid45 = "60045";
00087     const char* ETSolid92 = "60042";
00088     const char* ETSolid95 = "60095";
00089 
00090     // set Material id # to be used as default for all elements
00091     // will need to be subsequently reassigned inside ANSYS
00092     // Can, although have not, declare similar defaults for other attributes
00093     const char* MATDefault = "1";
00094 
00095     // create file streams for writing
00096     std::ofstream node_file;
00097     std::ofstream elem_file;
00098     std::ofstream ans_file;
00099 
00100     // get base filename from filename.ans
00101     std::string temp_string;
00102     std::string base_string;
00103     base_string.assign( file_name );
00104     base_string.replace( base_string.find_last_of( ".ans" ) - 3, 4, "" );
00105 
00106     // open node file for writing
00107     temp_string = base_string + ".node";
00108     node_file.open( temp_string.c_str() );
00109     node_file.setf( std::ios::scientific, std::ios::floatfield );
00110     node_file.precision( 13 );
00111 
00112     // open elem file for writing
00113     temp_string = base_string + ".elem";
00114     elem_file.open( temp_string.c_str() );
00115 
00116     // open ans file for writing
00117     ans_file.open( file_name );
00118     ans_file << "/prep7" << std::endl;
00119 
00120     // gather single output set
00121     EntityHandle output_set = 0;
00122     if( ent_handles && num_sets > 0 )
00123     {
00124         for( int i = 0; i < num_sets; i++ )
00125         {
00126             // from template, maybe can be removed
00127             result = mbImpl->unite_meshset( output_set, ent_handles[i] );
00128             if( result != MB_SUCCESS ) return result;
00129         }
00130     }
00131 
00132     // search for all nodes
00133     Range node_range;
00134     result = mbImpl->get_entities_by_type( output_set, MBVERTEX, node_range, true );
00135     if( result != MB_SUCCESS ) return result;
00136 
00137     // Commented out until Seg Fault taken care of in gather_nodes...
00138     // get any missing nodes which are needed for elements
00139     // Range all_ent_range,missing_range;
00140     // result=mbImpl->get_entities_by_handle(output_set,all_ent_range,true);
00141     // if(result !=MB_SUCCESS) return result;
00142     // result=mWriteIface->gather_nodes_from_elements(all_ent_range,0,missing_range);
00143     // node_range.merge(missing_range);
00144 
00145     // write the nodes
00146     double coord[3];
00147     for( Range::iterator it = node_range.begin(); it != node_range.end(); ++it )
00148     {
00149         EntityHandle node_handle = *it;
00150 
00151         result = mbImpl->get_coords( &node_handle, 1, coord );
00152         if( result != MB_SUCCESS ) return result;
00153 
00154         node_file.width( 8 );
00155         node_file << mbImpl->id_from_handle( node_handle );
00156         node_file.width( 20 );
00157         node_file << coord[0];
00158         node_file.width( 20 );
00159         node_file << coord[1];
00160         node_file.width( 20 );
00161         node_file << coord[2] << std::endl;
00162     }
00163 
00164     // update header to load nodes
00165     ans_file << "nread," << base_string << ",node" << std::endl;
00166 
00167     // search for all node sets (Dirichlet Sets)
00168     Range node_mesh_sets;
00169     int ns_id;
00170     result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, node_mesh_sets );
00171     if( result != MB_SUCCESS ) return result;
00172 
00173     for( Range::iterator ns_it = node_mesh_sets.begin(); ns_it != node_mesh_sets.end(); ++ns_it )
00174     {
00175         result = mbImpl->tag_get_data( mDirichletSetTag, &( *ns_it ), 1, &ns_id );
00176         if( result != MB_SUCCESS ) return result;
00177         std::vector< EntityHandle > node_vector;
00178         result = mbImpl->get_entities_by_handle( *ns_it, node_vector, true );
00179         if( result != MB_SUCCESS ) return result;
00180         // for every nodeset found, cycle through nodes in set:
00181         for( std::vector< EntityHandle >::iterator node_it = node_vector.begin(); node_it != node_vector.end();
00182              ++node_it )
00183         {
00184             int ns_node_id = mbImpl->id_from_handle( *node_it );
00185             if( node_it == node_vector.begin() )
00186             {
00187                 // select first node in new list
00188                 ans_file << "nsel,s,node,," << std::setw( 8 ) << ns_node_id << std::endl;
00189             }
00190             else
00191             {
00192                 // append node to list
00193                 ans_file << "nsel,a,node,," << std::setw( 8 ) << ns_node_id << std::endl;
00194             }
00195         }
00196         // create NS(#) nodeset
00197         ans_file << "cm,NS" << ns_id << ",node" << std::endl;
00198     }
00199 
00200     // ANSYS Element format:
00201     // I, J, K, L, M, N, O, P,etc... MAT, TYPE, REAL, SECNUM, ESYS, IEL
00202     // I-P are nodes of element
00203     // MAT = material number
00204     // TYPE = Element type number
00205     // REAL = Real constant set number
00206     // SECNUM = section attribute number
00207     // ESYS = coordinate system for nodes
00208     // IEL = element # (unique?)
00209     // For all nodes past 8, write on second line
00210 
00211     // Write all MBTET elements
00212     Range tet_range;
00213     result = mbImpl->get_entities_by_type( output_set, MBTET, tet_range, true );
00214     if( result != MB_SUCCESS ) return result;
00215     for( Range::iterator elem_it = tet_range.begin(); elem_it != tet_range.end(); ++elem_it )
00216     {
00217         EntityHandle elem_handle = *elem_it;
00218         int elem_id              = mbImpl->id_from_handle( elem_handle );
00219         std::vector< EntityHandle > conn;
00220         result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
00221         if( result != MB_SUCCESS ) return result;
00222         // make sure 4 or 10 node tet
00223         if( conn.size() != 4 && conn.size() != 10 )
00224         {
00225             std::cout << "Support not added for element type. \n";
00226             return MB_FAILURE;
00227         }
00228         // write information for 4 node tet
00229         if( conn.size() == 4 )
00230         {
00231             elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
00232             elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[2];
00233             elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];
00234             elem_file << std::setw( 8 ) << conn[3] << std::setw( 8 ) << conn[3];
00235 
00236             elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
00237             elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
00238             elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
00239             elem_file << std::endl;
00240         }
00241 
00242         // write information for 10 node tet
00243         if( conn.size() == 10 )
00244         {
00245             elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
00246             elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
00247             elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
00248             elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];
00249 
00250             elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid92;
00251             elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
00252             elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
00253             elem_file << std::endl;
00254 
00255             elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[9];
00256             elem_file << std::endl;
00257         }
00258     }
00259 
00260     // Write all MBHEX elements
00261     Range hex_range;
00262     result = mbImpl->get_entities_by_type( output_set, MBHEX, hex_range, true );
00263     if( result != MB_SUCCESS ) return result;
00264     for( Range::iterator elem_it = hex_range.begin(); elem_it != hex_range.end(); ++elem_it )
00265     {
00266         EntityHandle elem_handle = *elem_it;
00267         int elem_id              = mbImpl->id_from_handle( elem_handle );
00268         std::vector< EntityHandle > conn;
00269         result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
00270         if( result != MB_SUCCESS ) return result;
00271         // make sure supported hex type
00272         if( conn.size() != 8 && conn.size() != 20 )
00273         {
00274             std::cout << "Support not added for element type. \n";
00275             return MB_FAILURE;
00276         }
00277 
00278         // write information for 8 node hex
00279         if( conn.size() == 8 )
00280         {
00281             elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[1];
00282             elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
00283             elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
00284             elem_file << std::setw( 8 ) << conn[6] << std::setw( 8 ) << conn[7];
00285 
00286             elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
00287             elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
00288             elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
00289             elem_file << std::endl;
00290         }
00291 
00292         // write information for 20 node hex
00293         if( conn.size() == 20 )
00294         {
00295 
00296             elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[5];
00297             elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[0];
00298             elem_file << std::setw( 8 ) << conn[7] << std::setw( 8 ) << conn[6];
00299             elem_file << std::setw( 8 ) << conn[2] << std::setw( 8 ) << conn[3];
00300 
00301             elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid95;
00302             elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
00303             elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
00304             elem_file << std::endl;
00305 
00306             elem_file << std::setw( 8 ) << conn[16] << std::setw( 8 ) << conn[13];
00307             elem_file << std::setw( 8 ) << conn[8] << std::setw( 8 ) << conn[12];
00308             elem_file << std::setw( 8 ) << conn[18] << std::setw( 8 ) << conn[14];
00309             elem_file << std::setw( 8 ) << conn[10] << std::setw( 8 ) << conn[15];
00310             elem_file << std::setw( 8 ) << conn[19] << std::setw( 8 ) << conn[17];
00311             elem_file << std::setw( 8 ) << conn[9] << std::setw( 8 ) << conn[11];
00312             elem_file << std::endl;
00313         }
00314     }
00315     // Write all MBPRISM elements
00316     Range prism_range;
00317     result = mbImpl->get_entities_by_type( output_set, MBPRISM, prism_range, true );
00318     if( result != MB_SUCCESS ) return result;
00319     for( Range::iterator elem_it = prism_range.begin(); elem_it != prism_range.end(); ++elem_it )
00320     {
00321         EntityHandle elem_handle = *elem_it;
00322         int elem_id              = mbImpl->id_from_handle( elem_handle );
00323         std::vector< EntityHandle > conn;
00324         result = mbImpl->get_connectivity( &elem_handle, 1, conn, false );
00325         if( result != MB_SUCCESS ) return result;
00326         // make sure supported prism type
00327         if( conn.size() != 6 )
00328         {
00329             std::cout << "Support not added for element type. \n";
00330             return MB_FAILURE;
00331         }
00332 
00333         // write information for 6 node prism
00334         if( conn.size() == 6 )
00335         {
00336             elem_file << std::setw( 8 ) << conn[0] << std::setw( 8 ) << conn[3];
00337             elem_file << std::setw( 8 ) << conn[4] << std::setw( 8 ) << conn[4];
00338             elem_file << std::setw( 8 ) << conn[1] << std::setw( 8 ) << conn[2];
00339             elem_file << std::setw( 8 ) << conn[5] << std::setw( 8 ) << conn[5];
00340 
00341             elem_file << std::setw( 8 ) << MATDefault << std::setw( 8 ) << ETSolid45;
00342             elem_file << std::setw( 8 ) << "1" << std::setw( 8 ) << "1";
00343             elem_file << std::setw( 8 ) << "0" << std::setw( 8 ) << elem_id;
00344             elem_file << std::endl;
00345         }
00346     }
00347 
00348     // create element types (for now writes all, even if not used)
00349     ans_file << "et," << ETSolid45 << ",SOLID45" << std::endl;
00350     ans_file << "et," << ETSolid92 << ",SOLID92" << std::endl;
00351     ans_file << "et," << ETSolid95 << ",SOLID95" << std::endl;
00352 
00353     // xxx pyramids, other elements later...
00354 
00355     // write header to load elements
00356     ans_file << "eread," << base_string << ",elem" << std::endl;
00357 
00358     // search for all side sets (Neumann)
00359     Range side_mesh_sets;
00360     int ss_id;
00361     result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, side_mesh_sets );
00362     if( result != MB_SUCCESS ) return result;
00363     // cycle through all sets found
00364     for( Range::iterator ss_it = side_mesh_sets.begin(); ss_it != side_mesh_sets.end(); ++ss_it )
00365     {
00366         result = mbImpl->tag_get_data( mNeumannSetTag, &( *ss_it ), 1, &ss_id );
00367         if( result != MB_SUCCESS ) return result;
00368         std::vector< EntityHandle > elem_vector;
00369         result = mbImpl->get_entities_by_handle( *ss_it, elem_vector, true );
00370         if( result != MB_SUCCESS ) return result;
00371 
00372         // cycle through elements in current side set
00373         for( std::vector< EntityHandle >::iterator elem_it = elem_vector.begin(); elem_it != elem_vector.end();
00374              ++elem_it )
00375         {
00376             EntityHandle elem_handle = *elem_it;
00377 
00378             // instead of selecting current element in set, select its nodes...
00379             std::vector< EntityHandle > conn;
00380             result = mbImpl->get_connectivity( &elem_handle, 1, conn );
00381             if( result != MB_SUCCESS ) return result;
00382             if( elem_it == elem_vector.begin() )
00383             {
00384                 ans_file << "nsel,s,node,," << std::setw( 8 ) << conn[0] << std::endl;
00385                 for( unsigned int i = 1; i < conn.size(); i++ )
00386                 {
00387                     ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
00388                 }
00389             }
00390             else
00391             {
00392                 for( unsigned int i = 0; i < conn.size(); i++ )
00393                 {
00394                     ans_file << "nsel,a,node,," << std::setw( 8 ) << conn[i] << std::endl;
00395                 }
00396             }
00397         }
00398         // create SS(#) node set
00399         ans_file << "cm,SS" << ss_id << ",node" << std::endl;
00400     }
00401 
00402     // Gather all element blocks
00403     Range matset;
00404     int mat_id;
00405     result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, matset );
00406     if( result != MB_SUCCESS ) return result;
00407     // cycle through all elem blocks
00408     for( Range::iterator mat_it = matset.begin(); mat_it != matset.end(); ++mat_it )
00409     {
00410         EntityHandle matset_handle = *mat_it;
00411         result                     = mbImpl->tag_get_data( mMaterialSetTag, &matset_handle, 1, &mat_id );
00412         if( result != MB_SUCCESS ) return result;
00413         std::vector< EntityHandle > mat_vector;
00414         result = mbImpl->get_entities_by_handle( *mat_it, mat_vector, true );
00415         if( result != MB_SUCCESS ) return result;
00416         // cycle through elements in current mat set
00417         for( std::vector< EntityHandle >::iterator elem_it = mat_vector.begin(); elem_it != mat_vector.end();
00418              ++elem_it )
00419         {
00420             EntityHandle elem_handle = *elem_it;
00421             int elem_id              = mbImpl->id_from_handle( elem_handle );
00422             if( elem_it == mat_vector.begin() )
00423             { ans_file << "esel,s,elem,," << std::setw( 8 ) << elem_id << std::endl; }
00424             else
00425             {
00426                 ans_file << "esel,a,elem,," << std::setw( 8 ) << elem_id << std::endl;
00427             }
00428         }
00429         // for each matset, write block command
00430         ans_file << "cm,EB" << mat_id << ",elem" << std::endl;
00431     }
00432 
00433     // close all file streams
00434     node_file.close();
00435     elem_file.close();
00436     ans_file.close();
00437 
00438     return MB_SUCCESS;
00439 }
00440 
00441 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines