1
15
16 #ifdef WIN32
17 #ifdef _DEBUG
18
19
20 #pragma warning( disable : 4786 )
21 #endif
22 #endif
23
24 #include "WriteVtk.hpp"
25 #include "moab/VtkUtil.hpp"
26 #include "SysUtil.hpp"
27
28 #include <fstream>
29 #include <iostream>
30 #include <cstdio>
31 #include <cassert>
32 #include <vector>
33 #include <set>
34 #include <map>
35 #include <iterator>
36
37 #include "moab/Interface.hpp"
38 #include "moab/Range.hpp"
39 #include "moab/CN.hpp"
40 #include "MBTagConventions.hpp"
41 #include "moab/WriteUtilIface.hpp"
42 #include "Internals.hpp"
43 #include "moab/FileOptions.hpp"
44
45 namespace moab
46 {
47
48 const int DEFAULT_PRECISION = 10;
49 const bool DEFAULT_STRICT = true;
50
51 WriterIface* WriteVtk::factory( Interface* iface )
52 {
53 return new WriteVtk( iface );
54 }
55
56 WriteVtk::WriteVtk( Interface* impl )
57 : mbImpl( impl ), writeTool( 0 ), mStrict( DEFAULT_STRICT ), freeNodes( 0 ), createOneNodeCells( false )
58 {
59 assert( impl != NULL );
60 impl->query_interface( writeTool );
61 }
62
63 WriteVtk::~WriteVtk()
64 {
65 mbImpl->release_interface( writeTool );
66 }
67
68 ErrorCode WriteVtk::write_file( const char* file_name,
69 const bool overwrite,
70 const FileOptions& opts,
71 const EntityHandle* output_list,
72 const int num_sets,
73 const std::vector< std::string >& ,
74 const Tag* tag_list,
75 int num_tags,
76 int )
77 {
78 ErrorCode rval;
79
80
81 int precision;
82 if( MB_SUCCESS != opts.get_int_option( "PRECISION", precision ) ) precision = DEFAULT_PRECISION;
83
84 if( MB_SUCCESS == opts.get_null_option( "STRICT" ) )
85 mStrict = true;
86 else if( MB_SUCCESS == opts.get_null_option( "RELAXED" ) )
87 mStrict = false;
88 else
89 mStrict = DEFAULT_STRICT;
90
91 if( MB_SUCCESS == opts.get_null_option( "CREATE_ONE_NODE_CELLS" ) ) createOneNodeCells = true;
92
93
94 Range nodes, elems;
95 rval = gather_mesh( output_list, num_sets, nodes, elems );
96 if( MB_SUCCESS != rval ) return rval;
97
98
99 if( !overwrite )
100 {
101 rval = writeTool->check_doesnt_exist( file_name );
102 if( MB_SUCCESS != rval ) return rval;
103 }
104
105
106 std::ofstream file( file_name );
107 if( !file )
108 {
109 MB_SET_ERR( MB_FILE_WRITE_ERROR, "Could not open file: " << file_name );
110 }
111 file.precision( precision );
112
113
114 if( ( rval = write_header( file ) ) != MB_SUCCESS || ( rval = write_nodes( file, nodes ) ) != MB_SUCCESS ||
115 ( rval = write_elems( file, nodes, elems ) ) != MB_SUCCESS ||
116 ( rval = write_tags( file, true, nodes, tag_list, num_tags ) ) != MB_SUCCESS ||
117 ( rval = write_tags( file, false, elems, tag_list, num_tags ) ) != MB_SUCCESS )
118 {
119 file.close();
120 remove( file_name );
121 return rval;
122 }
123
124 return MB_SUCCESS;
125 }
126
127 ErrorCode WriteVtk::gather_mesh( const EntityHandle* set_list, int num_sets, Range& nodes, Range& elems )
128 {
129 ErrorCode rval;
130 int e;
131
132 if( !set_list || !num_sets )
133 {
134 Range a;
135 rval = mbImpl->get_entities_by_handle( 0, a );
136 if( MB_SUCCESS != rval ) return rval;
137
138 Range::const_iterator node_i, elem_i, set_i;
139 node_i = a.lower_bound( a.begin(), a.end(), CREATE_HANDLE( MBVERTEX, 0, e ) );
140 elem_i = a.lower_bound( node_i, a.end(), CREATE_HANDLE( MBEDGE, 0, e ) );
141 set_i = a.lower_bound( elem_i, a.end(), CREATE_HANDLE( MBENTITYSET, 0, e ) );
142 nodes.merge( node_i, elem_i );
143 elems.merge( elem_i, set_i );
144
145
146 EntityType et = MBEDGE;
147 for( et++; et < MBENTITYSET; et++ )
148 {
149 if( VtkUtil::get_vtk_type( et, CN::VerticesPerEntity( et ) ) ) continue;
150 Range::iterator eit = elems.lower_bound( elems.begin(), elems.end(), CREATE_HANDLE( et, 0, e ) ),
151 ep1it = elems.lower_bound( elems.begin(), elems.end(), CREATE_HANDLE( et + 1, 0, e ) );
152 elems.erase( eit, ep1it );
153 }
154 }
155 else
156 {
157 std::set< EntityHandle > visited;
158 std::vector< EntityHandle > sets;
159 sets.reserve( num_sets );
160 std::copy( set_list, set_list + num_sets, std::back_inserter( sets ) );
161 while( !sets.empty() )
162 {
163
164 EntityHandle set = sets.back();
165 sets.pop_back();
166
167 if( !visited.insert( set ).second ) continue;
168
169 Range a;
170 rval = mbImpl->get_entities_by_handle( set, a );
171 if( MB_SUCCESS != rval ) return rval;
172
173 Range::const_iterator node_i, elem_i, set_i;
174 node_i = a.lower_bound( a.begin(), a.end(), CREATE_HANDLE( MBVERTEX, 0, e ) );
175 elem_i = a.lower_bound( node_i, a.end(), CREATE_HANDLE( MBEDGE, 0, e ) );
176 set_i = a.lower_bound( elem_i, a.end(), CREATE_HANDLE( MBENTITYSET, 0, e ) );
177 nodes.merge( node_i, elem_i );
178 elems.merge( elem_i, set_i );
179 std::copy( set_i, a.end(), std::back_inserter( sets ) );
180
181 a.clear();
182 rval = mbImpl->get_child_meshsets( set, a );
183 std::copy( a.begin(), a.end(), std::back_inserter( sets ) );
184 }
185
186 for( Range::const_iterator ei = elems.begin(); ei != elems.end(); ++ei )
187 {
188 std::vector< EntityHandle > connect;
189 rval = mbImpl->get_connectivity( &( *ei ), 1, connect );
190 if( MB_SUCCESS != rval ) return rval;
191
192 for( unsigned int i = 0; i < connect.size(); ++i )
193 nodes.insert( connect[i] );
194 }
195 }
196
197 if( nodes.empty() )
198 {
199 MB_SET_ERR( MB_ENTITY_NOT_FOUND, "Nothing to write" );
200 }
201
202 return MB_SUCCESS;
203 }
204
205 ErrorCode WriteVtk::write_header( std::ostream& stream )
206 {
207 stream << "# vtk DataFile Version 3.0" << std::endl;
208 stream << MOAB_VERSION_STRING << std::endl;
209 stream << "ASCII" << std::endl;
210 stream << "DATASET UNSTRUCTURED_GRID" << std::endl;
211 return MB_SUCCESS;
212 }
213
214 ErrorCode WriteVtk::write_nodes( std::ostream& stream, const Range& nodes )
215 {
216 ErrorCode rval;
217
218 stream << "POINTS " << nodes.size() << " double" << std::endl;
219
220 double coords[3];
221 for( Range::const_iterator i = nodes.begin(); i != nodes.end(); ++i )
222 {
223 coords[1] = coords[2] = 0.0;
224 rval = mbImpl->get_coords( &( *i ), 1, coords );
225 if( MB_SUCCESS != rval ) return rval;
226 stream << coords[0] << ' ' << coords[1] << ' ' << coords[2] << std::endl;
227 }
228
229 return MB_SUCCESS;
230 }
231
232 ErrorCode WriteVtk::write_elems( std::ostream& stream, const Range& nodes, const Range& elems )
233 {
234 ErrorCode rval;
235
236 Range connectivity;
237 rval = mbImpl->get_connectivity( elems, connectivity );MB_CHK_ERR( rval );
238
239 Range nodes_from_connectivity = connectivity.subset_by_type( MBVERTEX );
240 Range faces_from_connectivity =
241 subtract( connectivity, nodes_from_connectivity );
242
243 Range connected_nodes;
244 rval = mbImpl->get_connectivity( faces_from_connectivity, connected_nodes );MB_CHK_ERR( rval );
245 connected_nodes.merge( nodes_from_connectivity );
246
247 Range free_nodes = subtract( nodes, connected_nodes );
248
249
250 unsigned long num_elems, num_uses;
251 num_elems = num_uses = elems.size();
252
253 std::map< EntityHandle, int > sizeFieldsPolyhedra;
254
255 for( Range::const_iterator i = elems.begin(); i != elems.end(); ++i )
256 {
257 EntityType type = mbImpl->type_from_handle( *i );
258 if( !VtkUtil::get_vtk_type( type, CN::VerticesPerEntity( type ) ) ) continue;
259
260 EntityHandle elem = *i;
261 const EntityHandle* connect = NULL;
262 int conn_len = 0;
263
264 std::vector< EntityHandle > storage;
265 rval = mbImpl->get_connectivity( elem, connect, conn_len, false, &storage );MB_CHK_ERR( rval );
266
267 num_uses += conn_len;
268
269 if( TYPE_FROM_HANDLE( elem ) == MBPOLYHEDRON )
270 {
271 int numFields = 1;
272 for( int j = 0; j < conn_len; j++ )
273 {
274 const EntityHandle* conn = NULL;
275 int num_nd = 0;
276 rval = mbImpl->get_connectivity( connect[j], conn, num_nd );MB_CHK_ERR( rval );
277 numFields += num_nd + 1;
278 }
279 sizeFieldsPolyhedra[elem] = numFields;
280 num_uses += ( numFields - conn_len );
281 }
282 }
283 freeNodes = (int)free_nodes.size();
284 if( !createOneNodeCells ) freeNodes = 0;
285 stream << "CELLS " << num_elems + freeNodes << ' ' << num_uses + 2 * freeNodes << std::endl;
286
287
288 std::vector< int > conn_data;
289 std::vector< unsigned > vtk_types( elems.size() + freeNodes );
290 std::vector< unsigned >::iterator t = vtk_types.begin();
291 for( Range::const_iterator i = elems.begin(); i != elems.end(); ++i )
292 {
293
294 EntityHandle elem = *i;
295 EntityType type = TYPE_FROM_HANDLE( elem );
296
297
298 const EntityHandle* connect = NULL;
299 int conn_len = 0;
300
301 std::vector< EntityHandle > storage;
302 rval = mbImpl->get_connectivity( elem, connect, conn_len, false, &storage );MB_CHK_ERR( rval );
303
304
305 const VtkElemType* vtk_type = VtkUtil::get_vtk_type( type, conn_len );
306 if( !vtk_type )
307 {
308
309 vtk_type = VtkUtil::get_vtk_type( type, conn_len - 1 );
310 if( vtk_type )
311 conn_len--;
312 else
313 {
314 MB_SET_ERR( MB_FAILURE, "Vtk file format does not support elements of type "
315 << CN::EntityTypeName( type ) << " (" << (int)type << ") with " << conn_len
316 << " nodes" );
317 }
318 }
319
320
321 *t = vtk_type->vtk_type;
322 ++t;
323
324 if( type != MBPOLYHEDRON )
325 {
326
327 assert( conn_len > 0 );
328 conn_data.resize( conn_len );
329 for( int j = 0; j < conn_len; ++j )
330 conn_data[j] = nodes.index( connect[j] );
331
332
333 stream << conn_len;
334 if( vtk_type->node_order )
335 for( int k = 0; k < conn_len; ++k )
336 stream << ' ' << conn_data[vtk_type->node_order[k]];
337 else
338 for( int k = 0; k < conn_len; ++k )
339 stream << ' ' << conn_data[k];
340 stream << std::endl;
341 }
342 else
343 {
344
345 stream << sizeFieldsPolyhedra[elem] << " " << conn_len;
346 for( int k = 0; k < conn_len; k++ )
347 {
348 EntityHandle face = connect[k];
349 const EntityHandle* conn = NULL;
350 int num_nodes = 0;
351 rval = mbImpl->get_connectivity( face, conn, num_nodes );MB_CHK_ERR( rval );
352
353 conn_data.resize( num_nodes );
354 for( int j = 0; j < num_nodes; ++j )
355 conn_data[j] = nodes.index( conn[j] );
356
357 stream << ' ' << num_nodes;
358
359 for( int j = 0; j < num_nodes; ++j )
360 stream << ' ' << conn_data[j];
361 }
362 stream << std::endl;
363 }
364 }
365
366 if( createOneNodeCells )
367 for( Range::const_iterator v = free_nodes.begin(); v != free_nodes.end(); ++v, ++t )
368 {
369 EntityHandle node = *v;
370 stream << "1 " << nodes.index( node ) << std::endl;
371 *t = 1;
372 }
373
374
375 stream << "CELL_TYPES " << vtk_types.size() << std::endl;
376 for( std::vector< unsigned >::const_iterator i = vtk_types.begin(); i != vtk_types.end(); ++i )
377 stream << *i << std::endl;
378
379 return MB_SUCCESS;
380 }
381
382 ErrorCode WriteVtk::write_tags( std::ostream& stream,
383 bool nodes,
384 const Range& entities,
385 const Tag* tag_list,
386 int num_tags )
387 {
388 ErrorCode rval;
389
390
391
392
393
394
395
396
397 EntityType low_type, high_type;
398 if( nodes )
399 {
400 low_type = MBVERTEX;
401 high_type = MBEDGE;
402 }
403 else
404 {
405 low_type = MBEDGE;
406 high_type = MBENTITYSET;
407 }
408
409
410 std::vector< Tag > tags;
411 std::vector< Tag >::iterator i;
412 rval = writeTool->get_tag_list( tags, tag_list, num_tags, false );
413 if( MB_SUCCESS != rval ) return rval;
414
415
416 bool entities_have_tags = false;
417 for( i = tags.begin(); i != tags.end(); ++i )
418 {
419
420 DataType dtype;
421 rval = mbImpl->tag_get_data_type( *i, dtype );
422 if( MB_SUCCESS != rval ) return rval;
423 if( dtype == MB_TYPE_HANDLE ) continue;
424
425
426
427 if( mStrict )
428 {
429 int count;
430 rval = mbImpl->tag_get_length( *i, count );
431 if( MB_SUCCESS != rval ) return rval;
432 if( count < 1 || ( count > 4 && count != 9 ) ) continue;
433 }
434
435
436 Range tagged;
437 for( EntityType type = low_type; type < high_type; ++type )
438 {
439 Range tmp_tagged;
440 rval = mbImpl->get_entities_by_type_and_tag( 0, type, &( *i ), 0, 1, tmp_tagged );
441 if( MB_SUCCESS != rval ) return rval;
442 tmp_tagged = intersect( tmp_tagged, entities );
443 tagged.merge( tmp_tagged );
444 }
445
446
447 if( !tagged.empty() )
448 {
449
450
451
452 if( !entities_have_tags )
453 {
454 entities_have_tags = true;
455 if( nodes )
456 stream << "POINT_DATA " << entities.size() << std::endl;
457 else
458 stream << "CELL_DATA " << entities.size() + freeNodes << std::endl;
459 }
460
461
462 rval = write_tag( stream, *i, entities, tagged );
463 if( MB_SUCCESS != rval ) return rval;
464 }
465 }
466
467 return MB_SUCCESS;
468 }
469
470 template < typename T >
471 void WriteVtk::write_data( std::ostream& stream, const std::vector< T >& data, unsigned vals_per_tag )
472 {
473 typename std::vector< T >::const_iterator d = data.begin();
474 const unsigned long n = data.size() / vals_per_tag;
475
476 for( unsigned long i = 0; i < n; ++i )
477 {
478 for( unsigned j = 0; j < vals_per_tag; ++j, ++d )
479 {
480 if( sizeof( T ) == 1 )
481 stream << (unsigned int)*d << ' ';
482 else
483 stream << *d << ' ';
484 }
485 stream << std::endl;
486 }
487 }
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504 template < typename T >
505 ErrorCode WriteVtk::write_tag( std::ostream& stream, Tag tag, const Range& entities, const Range& tagged, const int )
506 {
507 ErrorCode rval;
508 int addFreeNodes = 0;
509 if( TYPE_FROM_HANDLE( entities[0] ) > MBVERTEX ) addFreeNodes = freeNodes;
510
511
512
513 const unsigned long n = entities.size() + addFreeNodes;
514
515
516
517 std::string name;
518 int vals_per_tag;
519 if( MB_SUCCESS != mbImpl->tag_get_name( tag, name ) || MB_SUCCESS != mbImpl->tag_get_length( tag, vals_per_tag ) )
520 return MB_FAILURE;
521
522
523
524
525 std::vector< T > data;
526 data.resize( n * vals_per_tag, 0 );
527
528 std::vector< T > def_value( vals_per_tag );
529 rval = mbImpl->tag_get_default_value( tag, &( def_value[0] ) );
530 if( MB_SUCCESS == rval ) SysUtil::setmem( &( data[0] ), &( def_value[0] ), vals_per_tag * sizeof( T ), n );
531
532 Range::const_iterator t = tagged.begin();
533 typename std::vector< T >::iterator d = data.begin();
534 for( Range::const_iterator i = entities.begin(); i != entities.end() && t != tagged.end(); ++i, d += vals_per_tag )
535 {
536 if( *i == *t )
537 {
538 ++t;
539 rval = mbImpl->tag_get_data( tag, &( *i ), 1, &( *d ) );
540 if( MB_SUCCESS != rval ) return rval;
541 }
542 }
543
544
545 write_data( stream, data, vals_per_tag );
546
547 return MB_SUCCESS;
548 }
549
550 ErrorCode WriteVtk::write_bit_tag( std::ostream& stream, Tag tag, const Range& entities, const Range& tagged )
551 {
552 ErrorCode rval;
553 const unsigned long n = entities.size();
554
555
556
557 std::string name;
558 int vals_per_tag;
559 if( MB_SUCCESS != mbImpl->tag_get_name( tag, name ) || MB_SUCCESS != mbImpl->tag_get_length( tag, vals_per_tag ) )
560 return MB_FAILURE;
561
562 if( vals_per_tag > 8 )
563 {
564 MB_SET_ERR( MB_FAILURE, "Invalid tag size for bit tag \"" << name << "\"" );
565 }
566
567
568
569
570
571
572 std::vector< unsigned short > data;
573 data.resize( n * vals_per_tag, 0 );
574 Range::const_iterator t = tagged.begin();
575 std::vector< unsigned short >::iterator d = data.begin();
576 for( Range::const_iterator i = entities.begin(); i != entities.end() && t != tagged.end(); ++i )
577 {
578 if( *i == *t )
579 {
580 ++t;
581 unsigned char value;
582 rval = mbImpl->tag_get_data( tag, &( *i ), 1, &value );
583 for( int j = 0; j < vals_per_tag; ++j, ++d )
584 *d = (unsigned short)( value & ( 1 << j ) ? 1 : 0 );
585 if( MB_SUCCESS != rval ) return rval;
586 }
587 else
588 {
589
590 d += vals_per_tag;
591 }
592 }
593
594
595 write_data( stream, data, vals_per_tag );
596
597 return MB_SUCCESS;
598 }
599
600 ErrorCode WriteVtk::write_tag( std::ostream& s, Tag tag, const Range& entities, const Range& tagged )
601 {
602
603 std::string name;
604 DataType type;
605 int vals_per_tag;
606 if( MB_SUCCESS != mbImpl->tag_get_name( tag, name ) || MB_SUCCESS != mbImpl->tag_get_length( tag, vals_per_tag ) ||
607 MB_SUCCESS != mbImpl->tag_get_data_type( tag, type ) )
608 return MB_FAILURE;
609
610
611 if( MB_TYPE_HANDLE == type ) return MB_FAILURE;
612
613
614
615 for( std::string::iterator i = name.begin(); i != name.end(); ++i )
616 {
617 if( isspace( *i ) || iscntrl( *i ) ) *i = '_';
618 }
619
620
621 if( 3 == vals_per_tag && MB_TYPE_DOUBLE == type )
622 s << "VECTORS " << name << ' ' << VtkUtil::vtkTypeNames[type] << std::endl;
623 else if( 9 == vals_per_tag )
624 s << "TENSORS " << name << ' ' << VtkUtil::vtkTypeNames[type] << std::endl;
625 else
626 s << "SCALARS " << name << ' ' << VtkUtil::vtkTypeNames[type] << ' ' << vals_per_tag << std::endl
627 << "LOOKUP_TABLE default" << std::endl;
628
629
630 switch( type )
631 {
632 case MB_TYPE_OPAQUE:
633 return write_tag< unsigned char >( s, tag, entities, tagged, 0 );
634 case MB_TYPE_INTEGER:
635 return write_tag< int >( s, tag, entities, tagged, 0 );
636 case MB_TYPE_DOUBLE:
637 return write_tag< double >( s, tag, entities, tagged, 0 );
638 case MB_TYPE_BIT:
639 return write_bit_tag( s, tag, entities, tagged );
640 default:
641 return MB_FAILURE;
642 }
643 }
644
645 }