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