 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "WriteGmsh.hpp"
2 #include "moab/CN.hpp"
3 #include "MBTagConventions.hpp"
4 #include "MBParallelConventions.h"
5 #include "moab/Interface.hpp"
6 #include "moab/Range.hpp"
7 #include "moab/WriteUtilIface.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
20 WriterIface* WriteGmsh::factory( Interface* iface )
21 {
22 return new WriteGmsh( iface );
23 }
24
25 WriteGmsh::WriteGmsh( Interface* impl ) : mbImpl( impl )
26 {
27 impl->query_interface( mWriteIface );
28 }
29
30 WriteGmsh::~WriteGmsh()
31 {
32 mbImpl->release_interface( mWriteIface );
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();
72 mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag );
73 if( global_id ) mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
74 mbImpl->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, prtn_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
1.9.1.