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