WriteAns.cpp

Go to the documentation of this file.

/**
 * 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 )
 {
 std::cout << "Support not added for element type. \n";
 return MB_FAILURE;
 }
 
 // write information for 6 node prism
 if( conn.size() == 6 )
 {
 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