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