MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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