1
15
16 #ifdef WIN32
17 #ifdef _DEBUG
18
19
20 #pragma warning( disable : 4786 )
21 #endif
22 #endif
23
24 #include "WriteTemplate.hpp"
25
26 #include <utility>
27 #include <algorithm>
28 #include <ctime>
29 #include <string>
30 #include <vector>
31 #include <cstdio>
32 #include <cstring>
33 #include <iostream>
34
35 #include "moab/Interface.hpp"
36 #include "moab/Range.hpp"
37 #include "moab/CN.hpp"
38 #include <cassert>
39 #include "Internals.hpp"
40 #include "ExoIIUtil.hpp"
41 #include "MBTagConventions.hpp"
42 #include "moab/WriteUtilIface.hpp"
43
44 namespace moab
45 {
46
47 WriterIface* WriteTemplate::factory( Interface* iface )
48 {
49 return new WriteTemplate( iface );
50 }
51
52 WriteTemplate::WriteTemplate( Interface* impl ) : mbImpl( impl )
53 {
54 assert( impl != NULL );
55
56 impl->query_interface( mWriteIface );
57
58
59
60 int negone = -1;
61 impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
62 &negone );
63
64 impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
65 &negone );
66
67 impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
68 &negone );
69
70 mGlobalIdTag = impl->globalId_tag();
71
72 impl->tag_get_handle( "WriteTemplate element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
73 }
74
75 WriteTemplate::~WriteTemplate()
76 {
77 mbImpl->release_interface( mWriteIface );
78 mbImpl->tag_delete( mEntityMark );
79 }
80
81 void WriteTemplate::reset_matset( std::vector< WriteTemplate::MaterialSetData >& matset_info )
82 {
83 std::vector< WriteTemplate::MaterialSetData >::iterator iter;
84
85 for( iter = matset_info.begin(); iter != matset_info.end(); ++iter )
86 delete( *iter ).elements;
87 }
88
89 ErrorCode WriteTemplate::write_file( const char* file_name,
90 const bool ,
91 const FileOptions& ,
92 const EntityHandle* ent_handles,
93 const int num_sets,
94 const std::vector< std::string >& ,
95 const Tag* ,
96 int ,
97 int )
98 {
99 assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
100
101
102 if( NULL == strstr( file_name, ".template" ) ) return MB_FAILURE;
103
104 std::vector< EntityHandle > matsets, dirsets, neusets;
105
106 fileName = file_name;
107
108
109
110 if( num_sets == 0 )
111 {
112
113 Range this_range;
114 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
115 std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
116 this_range.clear();
117 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
118 std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
119 this_range.clear();
120 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
121 std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
122 }
123 else
124 {
125 int dummy;
126 for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
127 {
128 if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
129 matsets.push_back( *iter );
130 else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
131 dirsets.push_back( *iter );
132 else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
133 neusets.push_back( *iter );
134 }
135 }
136
137
138 if( matsets.empty() && dirsets.empty() && neusets.empty() ) return MB_FILE_WRITE_ERROR;
139
140 std::vector< WriteTemplate::MaterialSetData > matset_info;
141 std::vector< WriteTemplate::DirichletSetData > dirset_info;
142 std::vector< WriteTemplate::NeumannSetData > neuset_info;
143
144 MeshInfo mesh_info;
145
146 matset_info.clear();
147 if( gather_mesh_information( mesh_info, matset_info, neuset_info, dirset_info, matsets, neusets, dirsets ) !=
148 MB_SUCCESS )
149 {
150 reset_matset( matset_info );
151 return MB_FAILURE;
152 }
153
154
155 if( false )
156 {
157 reset_matset( matset_info );
158 return MB_FAILURE;
159 }
160
161 if( initialize_file( mesh_info ) != MB_SUCCESS )
162 {
163 reset_matset( matset_info );
164 return MB_FAILURE;
165 }
166
167 if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
168 {
169 reset_matset( matset_info );
170 return MB_FAILURE;
171 }
172
173 if( write_matsets( mesh_info, matset_info, neuset_info ) )
174 {
175 reset_matset( matset_info );
176 return MB_FAILURE;
177 }
178
179 return MB_SUCCESS;
180 }
181
182 ErrorCode WriteTemplate::gather_mesh_information( MeshInfo& mesh_info,
183 std::vector< WriteTemplate::MaterialSetData >& matset_info,
184 std::vector< WriteTemplate::NeumannSetData >& neuset_info,
185 std::vector< WriteTemplate::DirichletSetData >& dirset_info,
186 std::vector< EntityHandle >& matsets,
187 std::vector< EntityHandle >& neusets,
188 std::vector< EntityHandle >& dirsets )
189 {
190 std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
191
192 mesh_info.num_nodes = 0;
193 mesh_info.num_elements = 0;
194 mesh_info.num_matsets = 0;
195
196 int id = 0;
197
198 vector_iter = matsets.begin();
199 end_vector_iter = matsets.end();
200
201 mesh_info.num_matsets = matsets.size();
202
203 std::vector< EntityHandle > parent_meshsets;
204
205
206 mbImpl->tag_delete( mEntityMark );
207 mbImpl->tag_get_handle( "WriteTemplate element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
208
209 int highest_dimension_of_element_matsets = 0;
210
211 for( vector_iter = matsets.begin(); vector_iter != matsets.end(); ++vector_iter )
212 {
213 WriteTemplate::MaterialSetData matset_data;
214 matset_data.elements = new Range;
215
216
217 if( mbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
218
219
220 Range dummy_range;
221 mbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
222
223
224 Range::iterator entity_iter = dummy_range.end();
225 --entity_iter;
226 int this_dim = CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) );
227 entity_iter = dummy_range.begin();
228 while( entity_iter != dummy_range.end() && CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) ) != this_dim )
229 ++entity_iter;
230
231 if( entity_iter != dummy_range.end() )
232 std::copy( entity_iter, dummy_range.end(), range_inserter( *( matset_data.elements ) ) );
233
234 assert( matset_data.elements->begin() == matset_data.elements->end() ||
235 CN::Dimension( TYPE_FROM_HANDLE( *( matset_data.elements->begin() ) ) ) == this_dim );
236
237
238 if( mbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
239 {
240 MB_SET_ERR( MB_FAILURE, "Couldn't get matset id from a tag for an element matset" );
241 }
242
243 matset_data.id = id;
244 matset_data.number_attributes = 0;
245
246
247 Range::iterator elem_range_iter, end_elem_range_iter;
248 elem_range_iter = matset_data.elements->begin();
249 end_elem_range_iter = matset_data.elements->end();
250
251
252
253 EntityType entity_type = TYPE_FROM_HANDLE( *elem_range_iter );
254 --end_elem_range_iter;
255 if( entity_type != TYPE_FROM_HANDLE( *( end_elem_range_iter++ ) ) )
256 {
257 MB_SET_ERR( MB_FAILURE, "Entities in matset " << id << " not of common type" );
258 }
259
260 int dimension = CN::Dimension( entity_type );
261
262 if( dimension > highest_dimension_of_element_matsets ) highest_dimension_of_element_matsets = dimension;
263
264 matset_data.moab_type = mbImpl->type_from_handle( *( matset_data.elements->begin() ) );
265 if( MBMAXTYPE == matset_data.moab_type ) return MB_FAILURE;
266
267 std::vector< EntityHandle > tmp_conn;
268 mbImpl->get_connectivity( &( *( matset_data.elements->begin() ) ), 1, tmp_conn );
269 matset_data.element_type =
270 ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
271
272 if( matset_data.element_type == EXOII_MAX_ELEM_TYPE )
273 {
274 MB_SET_ERR( MB_FAILURE, "Element type in matset " << id << " didn't get set correctly" );
275 }
276
277 matset_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[matset_data.element_type];
278
279
280 matset_data.number_elements = matset_data.elements->size();
281
282
283 mesh_info.num_elements += matset_data.number_elements;
284
285
286 mWriteIface->gather_nodes_from_elements( *matset_data.elements, mEntityMark, mesh_info.nodes );
287
288 if( !neusets.empty() )
289 {
290
291 for( Range::iterator iter = matset_data.elements->begin(); iter != matset_data.elements->end(); ++iter )
292 {
293 unsigned char bit = 0x1;
294 mbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
295 }
296 }
297
298 matset_info.push_back( matset_data );
299 }
300
301
302 if( mesh_info.num_dim == 0 )
303 {
304
305 if( highest_dimension_of_element_matsets < 2 )
306 mesh_info.num_dim = 3;
307 else
308 mesh_info.num_dim = highest_dimension_of_element_matsets;
309 }
310
311 Range::iterator range_iter, end_range_iter;
312 range_iter = mesh_info.nodes.begin();
313 end_range_iter = mesh_info.nodes.end();
314
315 mesh_info.num_nodes = mesh_info.nodes.size();
316
317
318
319 vector_iter = dirsets.begin();
320 end_vector_iter = dirsets.end();
321
322 for( ; vector_iter != end_vector_iter; ++vector_iter )
323 {
324 WriteTemplate::DirichletSetData dirset_data;
325 dirset_data.id = 0;
326 dirset_data.number_nodes = 0;
327
328
329 if( mbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
330 {
331 MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for dirset " << id );
332 }
333
334 dirset_data.id = id;
335
336 std::vector< EntityHandle > node_vector;
337
338 if( mbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
339 {
340 MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in dirset " << id );
341 }
342
343 std::vector< EntityHandle >::iterator iter, end_iter;
344 iter = node_vector.begin();
345 end_iter = node_vector.end();
346
347 unsigned char node_marked = 0;
348 ErrorCode result;
349 for( ; iter != end_iter; ++iter )
350 {
351 if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
352 result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
353
354 if( 0x1 == node_marked ) dirset_data.nodes.push_back( *iter );
355 }
356
357 dirset_data.number_nodes = dirset_data.nodes.size();
358 dirset_info.push_back( dirset_data );
359 }
360
361
362 vector_iter = neusets.begin();
363 end_vector_iter = neusets.end();
364
365 for( ; vector_iter != end_vector_iter; ++vector_iter )
366 {
367 WriteTemplate::NeumannSetData neuset_data;
368
369
370 if( mbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
371
372 neuset_data.id = id;
373 neuset_data.mesh_set_handle = *vector_iter;
374
375
376
377 Range forward_elems, reverse_elems;
378 if( get_neuset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
379
380 ErrorCode result = get_valid_sides( forward_elems, 1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
381 result = get_valid_sides( reverse_elems, -1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
382
383 neuset_data.number_elements = neuset_data.elements.size();
384 neuset_info.push_back( neuset_data );
385 }
386
387 return MB_SUCCESS;
388 }
389
390 ErrorCode WriteTemplate::get_valid_sides( Range& elems, const int sense, WriteTemplate::NeumannSetData& neuset_data )
391 {
392
393
394 unsigned char element_marked = 0;
395 ErrorCode result;
396 for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
397 {
398
399 result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
400
401 if( 0x1 == element_marked )
402 {
403 neuset_data.elements.push_back( *iter );
404
405
406 neuset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
407 }
408 else
409 {
410 std::vector< EntityHandle > parents;
411 int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
412
413
414 if( mbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
415 {
416 MB_SET_ERR( MB_FAILURE, "Couldn't get adjacencies for neuset" );
417 }
418
419 if( !parents.empty() )
420 {
421
422 for( unsigned int k = 0; k < parents.size(); k++ )
423 {
424 result = mbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
425
426 int side_no, this_sense, this_offset;
427 if( 0x1 == element_marked &&
428 mbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
429 this_sense == sense )
430 {
431 neuset_data.elements.push_back( parents[k] );
432 neuset_data.side_numbers.push_back( side_no + 1 );
433 break;
434 }
435 }
436 }
437 else
438 {
439 MB_SET_ERR( MB_FAILURE, "No parent element exists for element in neuset " << neuset_data.id );
440 }
441 }
442 }
443
444 return MB_SUCCESS;
445 }
446
447 ErrorCode WriteTemplate::write_nodes( const int num_nodes, const Range& nodes, const int dimension )
448 {
449
450 ErrorCode result;
451 Tag trans_tag;
452 result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
453 bool transform_needed = true;
454 if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
455
456 int num_coords_to_fill = transform_needed ? 3 : dimension;
457
458 std::vector< double* > coord_arrays( 3 );
459 coord_arrays[0] = new double[num_nodes];
460 coord_arrays[1] = new double[num_nodes];
461 coord_arrays[2] = NULL;
462
463 if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
464
465 result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 0, coord_arrays );
466 if( result != MB_SUCCESS )
467 {
468 delete[] coord_arrays[0];
469 delete[] coord_arrays[1];
470 if( coord_arrays[2] ) delete[] coord_arrays[2];
471 return result;
472 }
473
474 if( transform_needed )
475 {
476 double trans_matrix[16];
477 const EntityHandle mesh = 0;
478 result = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
479
480 for( int i = 0; i < num_nodes; i++ )
481 {
482 double vec1[3];
483 double vec2[3];
484
485 vec2[0] = coord_arrays[0][i];
486 vec2[1] = coord_arrays[1][i];
487 vec2[2] = coord_arrays[2][i];
488
489 for( int row = 0; row < 3; row++ )
490 {
491 vec1[row] = 0.0;
492 for( int col = 0; col < 3; col++ )
493 vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
494 }
495
496 coord_arrays[0][i] = vec1[0];
497 coord_arrays[1][i] = vec1[1];
498 coord_arrays[2][i] = vec1[2];
499 }
500 }
501
502
503
504
505
506
507 delete[] coord_arrays[0];
508 delete[] coord_arrays[1];
509 if( coord_arrays[2] ) delete[] coord_arrays[2];
510
511 return MB_SUCCESS;
512 }
513
514 ErrorCode WriteTemplate::write_matsets(
515 MeshInfo& ,
516 std::vector< WriteTemplate::MaterialSetData >& matset_data,
517 std::vector< WriteTemplate::NeumannSetData >& )
518 {
519 unsigned int i;
520 std::vector< int > connect;
521 const EntityHandle* connecth;
522 int num_connecth;
523 ErrorCode result;
524
525
526 connect.reserve( 31 );
527 Range::iterator rit;
528
529 WriteTemplate::MaterialSetData matset;
530 for( i = 0; i < matset_data.size(); i++ )
531 {
532 matset = matset_data[i];
533
534 for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
535 {
536
537 result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
538 if( MB_SUCCESS != result ) return result;
539
540
541 result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[0] );
542 if( MB_SUCCESS != result ) return result;
543
544
545
546
547 if( false ) return MB_FAILURE;
548 }
549 }
550
551 return MB_SUCCESS;
552 }
553
554 ErrorCode WriteTemplate::initialize_file( MeshInfo& mesh_info )
555 {
556
557
558 int coord_size, ncoords;
559
560 coord_size = mesh_info.num_dim;
561 std::cout << "Coord_size = " << coord_size << std::endl;
562
563
564 ncoords = mesh_info.num_nodes;
565 std::cout << "ncoords = " << ncoords << std::endl;
566
567
568
570
571
572
573
574 return MB_SUCCESS;
575 }
576
577 ErrorCode WriteTemplate::open_file( const char* filename )
578 {
579
580 if( strlen( (const char*)filename ) == 0 )
581 {
582 MB_SET_ERR( MB_FAILURE, "Output filename not specified" );
583 }
584
585
586
587
588 if( false )
589 {
590 MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
591 }
592
593 return MB_SUCCESS;
594 }
595
596 ErrorCode WriteTemplate::get_neuset_elems( EntityHandle neuset,
597 int current_sense,
598 Range& forward_elems,
599 Range& reverse_elems )
600 {
601 Range neuset_elems, neuset_meshsets;
602
603
604
605 Tag sense_tag = 0;
606 mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
607
608
609 ErrorCode result = mbImpl->get_entities_by_handle( neuset, neuset_elems, true );
610 if( MB_FAILURE == result ) return result;
611
612
613 Range::iterator range_iter = neuset_elems.begin();
614 while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != neuset_elems.end() )
615 ++range_iter;
616
617
618 if( range_iter != neuset_elems.end() )
619 {
620 std::copy( range_iter, neuset_elems.end(), range_inserter( neuset_meshsets ) );
621 neuset_elems.erase( range_iter, neuset_elems.end() );
622 }
623
624
625
626
627
628 Range::iterator dum_it = neuset_elems.end();
629 --dum_it;
630 int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
631 dum_it = neuset_elems.begin();
632 while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != neuset_elems.end() )
633 ++dum_it;
634
635 if( current_sense == 1 || current_sense == 0 )
636 std::copy( dum_it, neuset_elems.end(), range_inserter( forward_elems ) );
637 if( current_sense == -1 || current_sense == 0 )
638 std::copy( dum_it, neuset_elems.end(), range_inserter( reverse_elems ) );
639
640
641
642 for( range_iter = neuset_meshsets.begin(); range_iter != neuset_meshsets.end(); ++range_iter )
643 {
644
645 int this_sense;
646 if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
647 this_sense = 1;
648
649
650 get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
651 }
652
653 return result;
654 }
655
656 }