MOAB: Mesh Oriented datABase  (version 5.2.1)
WriteGmsh.cpp
Go to the documentation of this file.
00001 #include "WriteGmsh.hpp"
00002 #include "moab/CN.hpp"
00003 #include "MBTagConventions.hpp"
00004 #include "MBParallelConventions.h"
00005 #include "moab/Interface.hpp"
00006 #include "moab/Range.hpp"
00007 #include "moab/WriteUtilIface.hpp"
00008 #include "moab/FileOptions.hpp"
00009 #include "GmshUtil.hpp"
00010 
00011 #include <fstream>
00012 #include <map>
00013 #include <set>
00014 
00015 namespace moab
00016 {
00017 
00018 const int DEFAULT_PRECISION = 10;
00019 
00020 WriterIface* WriteGmsh::factory( Interface* iface )
00021 {
00022     return new WriteGmsh( iface );
00023 }
00024 
00025 WriteGmsh::WriteGmsh( Interface* impl ) : mbImpl( impl )
00026 {
00027     impl->query_interface( mWriteIface );
00028 }
00029 
00030 WriteGmsh::~WriteGmsh()
00031 {
00032     mbImpl->release_interface( mWriteIface );
00033 }
00034 
00035 // A structure to store per-element information.
00036 struct ElemInfo
00037 {
00038     void set( int st, int idt )
00039     {
00040         while( count < st )
00041             sets[count++] = 0;
00042         sets[count++] = idt;
00043     }
00044     int count;    // Number of valid entries in sets[]
00045     int sets[3];  // IDs of owning block, geom, and partition; respectively
00046     int id;       // Global ID of element
00047     int type;     // Gmsh element type
00048 };
00049 
00050 //! Writes out a file
00051 ErrorCode WriteGmsh::write_file( const char* file_name, const bool overwrite, const FileOptions& options,
00052                                  const EntityHandle* output_list, const int num_sets,
00053                                  const std::vector< std::string >& /* qa_list */, const Tag* /* tag_list */,
00054                                  int /* num_tags */, int /* export_dimension */ )
00055 {
00056     ErrorCode rval;
00057     Tag global_id = 0, block_tag = 0, geom_tag = 0, prtn_tag = 0;
00058 
00059     if( !overwrite )
00060     {
00061         rval = mWriteIface->check_doesnt_exist( file_name );
00062         if( MB_SUCCESS != rval ) return rval;
00063     }
00064 
00065     // Get tags
00066     global_id = mbImpl->globalId_tag();
00067     mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag );
00068     if( global_id ) mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
00069     mbImpl->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, prtn_tag );
00070 
00071     // Define arrays to hold entity sets of interest
00072     Range sets[3];
00073     Tag set_tags[] = { block_tag, geom_tag, prtn_tag };
00074     Tag set_ids[]  = { block_tag, 0 /*global_id*/, prtn_tag };
00075 
00076     // Get entities to write
00077     Range elements, nodes;
00078     if( !output_list )
00079     {
00080         rval = mbImpl->get_entities_by_dimension( 0, 0, nodes, false );
00081         if( MB_SUCCESS != rval ) return rval;
00082         for( int d = 1; d <= 3; ++d )
00083         {
00084             Range tmp_range;
00085             rval = mbImpl->get_entities_by_dimension( 0, d, tmp_range, false );
00086             if( MB_SUCCESS != rval ) return rval;
00087             elements.merge( tmp_range );
00088         }
00089 
00090         for( int s = 0; s < 3; ++s )
00091         {
00092             if( set_tags[s] )
00093             {
00094                 rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, set_tags + s, 0, 1, sets[s] );
00095                 if( MB_SUCCESS != rval ) return rval;
00096             }
00097         }
00098     }
00099     else
00100     {
00101         for( int i = 0; i < num_sets; ++i )
00102         {
00103             EntityHandle set = output_list[i];
00104             for( int d = 1; d < 3; ++d )
00105             {
00106                 Range tmp_range, tmp_nodes;
00107                 rval = mbImpl->get_entities_by_dimension( set, d, tmp_range, true );
00108                 if( rval != MB_SUCCESS ) return rval;
00109                 elements.merge( tmp_range );
00110                 rval = mbImpl->get_adjacencies( tmp_range, set, false, tmp_nodes );
00111                 if( rval != MB_SUCCESS ) return rval;
00112                 nodes.merge( tmp_nodes );
00113             }
00114 
00115             for( int s = 0; s < 3; ++s )
00116             {
00117                 if( set_tags[s] )
00118                 {
00119                     Range tmp_range;
00120                     rval = mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, set_tags + s, 0, 1, tmp_range );
00121                     if( MB_SUCCESS != rval ) return rval;
00122                     sets[s].merge( tmp_range );
00123                     int junk;
00124                     rval = mbImpl->tag_get_data( set_tags[s], &set, 1, &junk );
00125                     if( MB_SUCCESS == rval ) sets[s].insert( set );
00126                 }
00127             }
00128         }
00129     }
00130 
00131     if( elements.empty() ) { MB_SET_ERR( MB_ENTITY_NOT_FOUND, "Nothing to write" ); }
00132 
00133     // Get global IDs for all elements.
00134     // First try to get from tag.  If tag is not defined or not set
00135     // for all elements, use handle value instead.
00136     std::vector< int > global_id_array( elements.size() );
00137     std::vector< int >::iterator id_iter;
00138     if( !global_id || MB_SUCCESS != mbImpl->tag_get_data( global_id, elements, &global_id_array[0] ) )
00139     {
00140         id_iter = global_id_array.begin();
00141         for( Range::iterator i = elements.begin(); i != elements.end(); ++i, ++id_iter )
00142             *id_iter = mbImpl->id_from_handle( *i );
00143     }
00144 
00145     // Figure out the maximum ID value so we know where to start allocating
00146     // new IDs when we encounter ID conflicts.
00147     int max_id = 0;
00148     for( id_iter = global_id_array.begin(); id_iter != global_id_array.end(); ++id_iter )
00149         if( *id_iter > max_id ) max_id = *id_iter;
00150 
00151     // Initialize ElemInfo struct for each element
00152     std::map< EntityHandle, ElemInfo > elem_sets;  // Per-element info
00153     std::set< int > elem_global_ids;               // Temporary for finding duplicate IDs
00154     id_iter = global_id_array.begin();
00155     // Iterate backwards to give highest-dimension entities first dibs for
00156     // a conflicting ID.
00157     for( Range::reverse_iterator i = elements.rbegin(); i != elements.rend(); ++i )
00158     {
00159         int id = *id_iter;
00160         ++id_iter;
00161         if( !elem_global_ids.insert( id ).second ) id = ++max_id;
00162 
00163         ElemInfo& ei = elem_sets[*i];
00164         ei.count     = 0;
00165         ei.id        = id;
00166 
00167         EntityType type = mbImpl->type_from_handle( *i );
00168         int num_vtx;
00169         const EntityHandle* conn;
00170         rval = mbImpl->get_connectivity( *i, conn, num_vtx );
00171         if( MB_SUCCESS != rval ) return rval;
00172 
00173         ei.type = GmshUtil::get_gmsh_type( type, num_vtx );
00174         if( ei.type < 0 )
00175         {
00176             MB_SET_ERR( MB_FILE_WRITE_ERROR, "Gmem file format does not support element of type "
00177                                                  << CN::EntityTypeName( type ) << " with " << num_vtx << " vertices" );
00178         }
00179     }
00180     // Don't need these any more, free memory.
00181     elem_global_ids.clear();
00182     global_id_array.clear();
00183 
00184     // For each material set, geometry set, or partition; store
00185     // the ID of the set on each element.
00186     for( int s = 0; s < 3; ++s )
00187     {
00188         if( !set_tags[s] ) continue;
00189 
00190         for( Range::iterator i = sets[s].begin(); i != sets[s].end(); ++i )
00191         {
00192             int id;
00193             if( set_ids[s] )
00194             {
00195                 rval = mbImpl->tag_get_data( set_ids[s], &*i, 1, &id );
00196                 if( MB_SUCCESS != rval ) return rval;
00197             }
00198             else
00199                 id = mbImpl->id_from_handle( *i );
00200 
00201             Range elems;
00202             rval = mbImpl->get_entities_by_handle( *i, elems );
00203             if( MB_SUCCESS != rval ) return rval;
00204 
00205             elems = intersect( elems, elements );
00206             for( Range::iterator j = elems.begin(); j != elems.end(); ++j )
00207                 elem_sets[*j].set( s, id );
00208         }
00209     }
00210 
00211     // Create file
00212     std::ofstream out( file_name );
00213     if( !out ) return MB_FILE_DOES_NOT_EXIST;
00214 
00215     // Write header
00216     out << "$MeshFormat" << std::endl;
00217     out << "2.0 0 " << sizeof( double ) << std::endl;
00218     out << "$EndMeshFormat" << std::endl;
00219 
00220     // Set precision for node coordinates
00221     int precision;
00222     if( MB_SUCCESS != options.get_int_option( "PRECISION", precision ) ) precision = DEFAULT_PRECISION;
00223     const int old_precision = out.precision();
00224     out.precision( precision );
00225 
00226     // Write nodes
00227     out << "$Nodes" << std::endl;
00228     out << nodes.size() << std::endl;
00229     std::vector< double > coords( 3 * nodes.size() );
00230     rval = mbImpl->get_coords( nodes, &coords[0] );
00231     if( MB_SUCCESS != rval ) return rval;
00232     std::vector< double >::iterator c = coords.begin();
00233     for( Range::iterator i = nodes.begin(); i != nodes.end(); ++i )
00234     {
00235         out << mbImpl->id_from_handle( *i );
00236         out << " " << *c;
00237         ++c;
00238         out << " " << *c;
00239         ++c;
00240         out << " " << *c;
00241         ++c;
00242         out << std::endl;
00243     }
00244     out << "$EndNodes" << std::endl;
00245     coords.clear();
00246 
00247     // Restore stream state
00248     out.precision( old_precision );
00249 
00250     // Write elements
00251     out << "$Elements" << std::endl;
00252     out << elem_sets.size() << std::endl;
00253     for( std::map< EntityHandle, ElemInfo >::iterator i = elem_sets.begin(); i != elem_sets.end(); ++i )
00254     {
00255         int num_vtx;
00256         const EntityHandle* conn;
00257         rval = mbImpl->get_connectivity( i->first, conn, num_vtx );
00258         if( MB_SUCCESS != rval ) return rval;
00259         out << i->second.id << ' ' << i->second.type << ' ' << i->second.count;
00260         for( int j = 0; j < i->second.count; ++j )
00261             out << ' ' << i->second.sets[j];
00262 
00263         const int* order = GmshUtil::gmshElemTypes[i->second.type].node_order;
00264 
00265         // Need to re-order vertices
00266         if( order )
00267         {
00268             for( int j = 0; j < num_vtx; ++j )
00269                 out << ' ' << mbImpl->id_from_handle( conn[order[j]] );
00270         }
00271         else
00272         {
00273             for( int j = 0; j < num_vtx; ++j )
00274                 out << ' ' << mbImpl->id_from_handle( conn[j] );
00275         }
00276         out << std::endl;
00277     }
00278     out << "$EndElements" << std::endl;
00279 
00280     // Done
00281     return MB_SUCCESS;
00282 }
00283 
00284 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines