Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
WriteGmsh.cpp
Go to the documentation of this file.
1 #include "WriteGmsh.hpp"
2 #include "moab/CN.hpp"
3 #include "MBTagConventions.hpp"
5 #include "moab/Interface.hpp"
6 #include "moab/Range.hpp"
8 #include "moab/FileOptions.hpp"
9 #include "GmshUtil.hpp"
10 
11 #include <fstream>
12 #include <map>
13 #include <set>
14 
15 namespace moab
16 {
17 
18 const int DEFAULT_PRECISION = 10;
19 
21 {
22  return new WriteGmsh( iface );
23 }
24 
25 WriteGmsh::WriteGmsh( Interface* impl ) : mbImpl( impl )
26 {
28 }
29 
31 {
33 }
34 
35 // A structure to store per-element information.
36 struct ElemInfo
37 {
38  void set( int st, int idt )
39  {
40  while( count < st )
41  sets[count++] = 0;
42  sets[count++] = idt;
43  }
44  int count; // Number of valid entries in sets[]
45  int sets[3]; // IDs of owning block, geom, and partition; respectively
46  int id; // Global ID of element
47  int type; // Gmsh element type
48 };
49 
50 //! Writes out a file
51 ErrorCode WriteGmsh::write_file( const char* file_name,
52  const bool overwrite,
53  const FileOptions& options,
54  const EntityHandle* output_list,
55  const int num_sets,
56  const std::vector< std::string >& /* qa_list */,
57  const Tag* /* tag_list */,
58  int /* num_tags */,
59  int /* export_dimension */ )
60 {
61  ErrorCode rval;
62  Tag global_id = 0, block_tag = 0, geom_tag = 0, prtn_tag = 0;
63 
64  if( !overwrite )
65  {
66  rval = mWriteIface->check_doesnt_exist( file_name );
67  if( MB_SUCCESS != rval ) return rval;
68  }
69 
70  // Get tags
71  global_id = mbImpl->globalId_tag();
73  if( global_id ) mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
75 
76  // Define arrays to hold entity sets of interest
77  Range sets[3];
78  Tag set_tags[] = { block_tag, geom_tag, prtn_tag };
79  Tag set_ids[] = { block_tag, 0 /*global_id*/, prtn_tag };
80 
81  // Get entities to write
82  Range elements, nodes;
83  if( !output_list )
84  {
85  rval = mbImpl->get_entities_by_dimension( 0, 0, nodes, false );
86  if( MB_SUCCESS != rval ) return rval;
87  for( int d = 1; d <= 3; ++d )
88  {
89  Range tmp_range;
90  rval = mbImpl->get_entities_by_dimension( 0, d, tmp_range, false );
91  if( MB_SUCCESS != rval ) return rval;
92  elements.merge( tmp_range );
93  }
94 
95  for( int s = 0; s < 3; ++s )
96  {
97  if( set_tags[s] )
98  {
99  rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, set_tags + s, 0, 1, sets[s] );
100  if( MB_SUCCESS != rval ) return rval;
101  }
102  }
103  }
104  else
105  {
106  for( int i = 0; i < num_sets; ++i )
107  {
108  EntityHandle set = output_list[i];
109  for( int d = 1; d < 3; ++d )
110  {
111  Range tmp_range, tmp_nodes;
112  rval = mbImpl->get_entities_by_dimension( set, d, tmp_range, true );
113  if( rval != MB_SUCCESS ) return rval;
114  elements.merge( tmp_range );
115  rval = mbImpl->get_adjacencies( tmp_range, set, false, tmp_nodes );
116  if( rval != MB_SUCCESS ) return rval;
117  nodes.merge( tmp_nodes );
118  }
119 
120  for( int s = 0; s < 3; ++s )
121  {
122  if( set_tags[s] )
123  {
124  Range tmp_range;
125  rval = mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, set_tags + s, 0, 1, tmp_range );
126  if( MB_SUCCESS != rval ) return rval;
127  sets[s].merge( tmp_range );
128  int junk;
129  rval = mbImpl->tag_get_data( set_tags[s], &set, 1, &junk );
130  if( MB_SUCCESS == rval ) sets[s].insert( set );
131  }
132  }
133  }
134  }
135 
136  if( elements.empty() )
137  {
138  MB_SET_ERR( MB_ENTITY_NOT_FOUND, "Nothing to write" );
139  }
140 
141  // Get global IDs for all elements.
142  // First try to get from tag. If tag is not defined or not set
143  // for all elements, use handle value instead.
144  std::vector< int > global_id_array( elements.size() );
145  std::vector< int >::iterator id_iter;
146  if( !global_id || MB_SUCCESS != mbImpl->tag_get_data( global_id, elements, &global_id_array[0] ) )
147  {
148  id_iter = global_id_array.begin();
149  for( Range::iterator i = elements.begin(); i != elements.end(); ++i, ++id_iter )
150  *id_iter = mbImpl->id_from_handle( *i );
151  }
152 
153  // Figure out the maximum ID value so we know where to start allocating
154  // new IDs when we encounter ID conflicts.
155  int max_id = 0;
156  for( id_iter = global_id_array.begin(); id_iter != global_id_array.end(); ++id_iter )
157  if( *id_iter > max_id ) max_id = *id_iter;
158 
159  // Initialize ElemInfo struct for each element
160  std::map< EntityHandle, ElemInfo > elem_sets; // Per-element info
161  std::set< int > elem_global_ids; // Temporary for finding duplicate IDs
162  id_iter = global_id_array.begin();
163  // Iterate backwards to give highest-dimension entities first dibs for
164  // a conflicting ID.
165  for( Range::reverse_iterator i = elements.rbegin(); i != elements.rend(); ++i )
166  {
167  int id = *id_iter;
168  ++id_iter;
169  if( !elem_global_ids.insert( id ).second ) id = ++max_id;
170 
171  ElemInfo& ei = elem_sets[*i];
172  ei.count = 0;
173  ei.id = id;
174 
175  EntityType type = mbImpl->type_from_handle( *i );
176  int num_vtx;
177  const EntityHandle* conn;
178  rval = mbImpl->get_connectivity( *i, conn, num_vtx );
179  if( MB_SUCCESS != rval ) return rval;
180 
181  ei.type = GmshUtil::get_gmsh_type( type, num_vtx );
182  if( ei.type < 0 )
183  {
184  MB_SET_ERR( MB_FILE_WRITE_ERROR, "Gmem file format does not support element of type "
185  << CN::EntityTypeName( type ) << " with " << num_vtx << " vertices" );
186  }
187  }
188  // Don't need these any more, free memory.
189  elem_global_ids.clear();
190  global_id_array.clear();
191 
192  // For each material set, geometry set, or partition; store
193  // the ID of the set on each element.
194  for( int s = 0; s < 3; ++s )
195  {
196  if( !set_tags[s] ) continue;
197 
198  for( Range::iterator i = sets[s].begin(); i != sets[s].end(); ++i )
199  {
200  int id;
201  if( set_ids[s] )
202  {
203  rval = mbImpl->tag_get_data( set_ids[s], &*i, 1, &id );
204  if( MB_SUCCESS != rval ) return rval;
205  }
206  else
207  id = mbImpl->id_from_handle( *i );
208 
209  Range elems;
210  rval = mbImpl->get_entities_by_handle( *i, elems );
211  if( MB_SUCCESS != rval ) return rval;
212 
213  elems = intersect( elems, elements );
214  for( Range::iterator j = elems.begin(); j != elems.end(); ++j )
215  elem_sets[*j].set( s, id );
216  }
217  }
218 
219  // Create file
220  std::ofstream out( file_name );
221  if( !out ) return MB_FILE_DOES_NOT_EXIST;
222 
223  // Write header
224  out << "$MeshFormat" << std::endl;
225  out << "2.0 0 " << sizeof( double ) << std::endl;
226  out << "$EndMeshFormat" << std::endl;
227 
228  // Set precision for node coordinates
229  int precision;
230  if( MB_SUCCESS != options.get_int_option( "PRECISION", precision ) ) precision = DEFAULT_PRECISION;
231  const int old_precision = out.precision();
232  out.precision( precision );
233 
234  // Write nodes
235  out << "$Nodes" << std::endl;
236  out << nodes.size() << std::endl;
237  std::vector< double > coords( 3 * nodes.size() );
238  rval = mbImpl->get_coords( nodes, &coords[0] );
239  if( MB_SUCCESS != rval ) return rval;
240  std::vector< double >::iterator c = coords.begin();
241  for( Range::iterator i = nodes.begin(); i != nodes.end(); ++i )
242  {
243  out << mbImpl->id_from_handle( *i );
244  out << " " << *c;
245  ++c;
246  out << " " << *c;
247  ++c;
248  out << " " << *c;
249  ++c;
250  out << std::endl;
251  }
252  out << "$EndNodes" << std::endl;
253  coords.clear();
254 
255  // Restore stream state
256  out.precision( old_precision );
257 
258  // Write elements
259  out << "$Elements" << std::endl;
260  out << elem_sets.size() << std::endl;
261  for( std::map< EntityHandle, ElemInfo >::iterator i = elem_sets.begin(); i != elem_sets.end(); ++i )
262  {
263  int num_vtx;
264  const EntityHandle* conn;
265  rval = mbImpl->get_connectivity( i->first, conn, num_vtx );
266  if( MB_SUCCESS != rval ) return rval;
267  out << i->second.id << ' ' << i->second.type << ' ' << i->second.count;
268  for( int j = 0; j < i->second.count; ++j )
269  out << ' ' << i->second.sets[j];
270 
271  const int* order = GmshUtil::gmshElemTypes[i->second.type].node_order;
272 
273  // Need to re-order vertices
274  if( order )
275  {
276  for( int j = 0; j < num_vtx; ++j )
277  out << ' ' << mbImpl->id_from_handle( conn[order[j]] );
278  }
279  else
280  {
281  for( int j = 0; j < num_vtx; ++j )
282  out << ' ' << mbImpl->id_from_handle( conn[j] );
283  }
284  out << std::endl;
285  }
286  out << "$EndElements" << std::endl;
287 
288  // Done
289  return MB_SUCCESS;
290 }
291 
292 } // namespace moab