 |
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 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 #ifdef WIN32
17 #ifdef _DEBUG
18 // turn off warnings that say they debugging identifier has been truncated
19 // this warning comes up when using some STL containers
20 #pragma warning( disable : 4786 )
21 #endif
22 #endif
23
24 #include "WriteNCDF.hpp"
25
26 #include "netcdf.h"
27 #include <utility>
28 #include <algorithm>
29 #include <ctime>
30 #include <string>
31 #include <vector>
32 #include <cstdio>
33 #include <cstring>
34 #include <cassert>
35
36 #include "moab/Interface.hpp"
37 #include "moab/Range.hpp"
38 #include "moab/CN.hpp"
39 #include "moab/FileOptions.hpp"
40 #include "MBTagConventions.hpp"
41 #include "Internals.hpp"
42 #include "ExoIIUtil.hpp"
43 #include "moab/WriteUtilIface.hpp"
44 #include "exodus_order.h"
45
46 #ifndef MOAB_HAVE_NETCDF
47 #error Attempt to compile WriteNCDF with NetCDF support disabled
48 #endif
49
50 #define CHAR_STR_LEN 128
51
52 namespace moab
53 {
54
55 const int TIME_STR_LEN = 11;
56
57 #define INS_ID( stringvar, prefix, id ) snprintf( stringvar, CHAR_STR_LEN, prefix, id )
58
59 #define GET_DIM( ncdim, name, val ) \
60 { \
61 int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) ); \
62 if( NC_NOERR == gdfail ) \
63 { \
64 size_t tmp_val; \
65 gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
66 if( NC_NOERR != gdfail ) \
67 { \
68 MB_SET_ERR( MB_FAILURE, "WriteNCDF:: couldn't get dimension length" ); \
69 } \
70 else \
71 ( val ) = tmp_val; \
72 } \
73 else \
74 ( val ) = 0; \
75 }
76
77 #define GET_DIMB( ncdim, name, varname, id, val ) \
78 INS_ID( name, CHAR_STR_LEN, varname, id ); \
79 GET_DIM( ncdim, name, val );
80
81 #define GET_VAR( name, id, dims ) \
82 { \
83 ( id ) = -1; \
84 int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
85 if( NC_NOERR == gvfail ) \
86 { \
87 int ndims; \
88 gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
89 if( NC_NOERR == gvfail ) \
90 { \
91 ( dims ).resize( ndims ); \
92 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
93 } \
94 } \
95 }
96
97 WriterIface* WriteNCDF::factory( Interface* iface )
98 {
99 return new WriteNCDF( iface );
100 }
101
102 WriteNCDF::WriteNCDF( Interface* impl )
103 : mdbImpl( impl ), ncFile( 0 ), mCurrentMeshHandle( 0 ), mGeomDimensionTag( 0 ), repeat_face_blocks( 0 )
104 {
105 assert( impl != NULL );
106
107 impl->query_interface( mWriteIface );
108
109 // Initialize in case tag_get_handle fails below
110 //! Get and cache predefined tag handles
111 int negone = -1;
112 impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
113 &negone );
114
115 impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
116 &negone );
117
118 impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
119 &negone );
120
121 mGlobalIdTag = impl->globalId_tag();
122
123 int dum_val_array[] = { -1, -1, -1, -1 };
124 impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag, MB_TAG_SPARSE | MB_TAG_CREAT,
125 dum_val_array );
126
127 impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
128 MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
129
130 impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag, MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
131
132 impl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
133 }
134
135 WriteNCDF::~WriteNCDF()
136 {
137 mdbImpl->release_interface( mWriteIface );
138
139 mdbImpl->tag_delete( mEntityMark );
140
141 if( 0 != ncFile ) ncFile = 0;
142 }
143
144 void WriteNCDF::reset_block( std::vector< MaterialSetData >& block_info )
145 {
146 std::vector< MaterialSetData >::iterator iter;
147
148 for( iter = block_info.begin(); iter != block_info.end(); ++iter )
149 {
150 iter->elements.clear();
151 }
152 }
153
154 void WriteNCDF::time_and_date( char* time_string, char* date_string )
155 {
156 struct tm* local_time;
157 time_t calendar_time;
158
159 calendar_time = time( NULL );
160 local_time = localtime( &calendar_time );
161
162 assert( NULL != time_string && NULL != date_string );
163
164 strftime( time_string, TIME_STR_LEN, "%H:%M:%S", local_time );
165 strftime( date_string, TIME_STR_LEN, "%m/%d/%Y", local_time );
166
167 // Terminate with NULL character
168 time_string[10] = '\0';
169 date_string[10] = '\0';
170 }
171
172 ErrorCode WriteNCDF::write_file( const char* exodus_file_name,
173 const bool overwrite,
174 const FileOptions& opts,
175 const EntityHandle* ent_handles,
176 const int num_sets,
177 const std::vector< std::string >& qa_records,
178 const Tag*,
179 int,
180 int user_dimension )
181 {
182 assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
183
184 if( user_dimension == 0 ) mdbImpl->get_dimension( user_dimension );
185
186 if( opts.get_null_option( "REPEAT_FACE_BLOCKS" ) == MB_SUCCESS ) repeat_face_blocks = 1;
187
188 std::vector< EntityHandle > blocks, nodesets, sidesets, entities;
189
190 // Separate into blocks, nodesets, sidesets
191
192 if( num_sets == 0 )
193 {
194 // Default to all defined block, nodeset and sideset-type sets
195 Range this_range;
196 mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
197 std::copy( this_range.begin(), this_range.end(), std::back_inserter( blocks ) );
198 this_range.clear();
199 mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
200 std::copy( this_range.begin(), this_range.end(), std::back_inserter( nodesets ) );
201 this_range.clear();
202 mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
203 std::copy( this_range.begin(), this_range.end(), std::back_inserter( sidesets ) );
204
205 // If there is nothing to write, write everything as one block.
206 if( blocks.empty() && nodesets.empty() && sidesets.empty() )
207 {
208 this_range.clear();
209 for( int d = user_dimension; d > 0 && this_range.empty(); --d )
210 mdbImpl->get_entities_by_dimension( 0, d, this_range, false );
211 if( this_range.empty() ) return MB_FILE_WRITE_ERROR;
212
213 EntityHandle block_handle;
214 int block_id = 1;
215 mdbImpl->create_meshset( MESHSET_SET, block_handle );
216 mdbImpl->tag_set_data( mMaterialSetTag, &block_handle, 1, &block_id );
217 mdbImpl->add_entities( block_handle, this_range );
218 blocks.push_back( block_handle );
219 }
220 }
221 else
222 {
223 int dummy;
224 for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
225 {
226 if( MB_SUCCESS == mdbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
227 blocks.push_back( *iter );
228 else if( MB_SUCCESS == mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
229 nodesets.push_back( *iter );
230 else if( MB_SUCCESS == mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
231 sidesets.push_back( *iter );
232 }
233 }
234
235 // If there is nothing to write just return.
236 if( blocks.empty() && nodesets.empty() && sidesets.empty() ) return MB_FILE_WRITE_ERROR;
237
238 // Try to get mesh information
239 ExodusMeshInfo mesh_info;
240
241 std::vector< MaterialSetData > block_info;
242 std::vector< NeumannSetData > sideset_info;
243 std::vector< DirichletSetData > nodeset_info;
244
245 mesh_info.num_dim = user_dimension;
246
247 if( qa_records.empty() )
248 {
249 // qa records are empty - initialize some MB-standard ones
250 mesh_info.qaRecords.push_back( "MB" );
251 mesh_info.qaRecords.push_back( "0.99" );
252 char string1[80], string2[80];
253 time_and_date( string2, string1 );
254 mesh_info.qaRecords.push_back( string2 );
255 mesh_info.qaRecords.push_back( string1 );
256 }
257 else
258 {
259 // Constrained to multiples of 4 qa records
260 assert( qa_records.size() % 4 == 0 );
261
262 std::copy( qa_records.begin(), qa_records.end(), std::back_inserter( mesh_info.qaRecords ) );
263 }
264
265 block_info.clear();
266 if( gather_mesh_information( mesh_info, block_info, sideset_info, nodeset_info, blocks, sidesets, nodesets ) !=
267 MB_SUCCESS )
268 {
269 reset_block( block_info );
270 return MB_FAILURE;
271 }
272
273 // Try to open the file after gather mesh info succeeds
274 int fail = nc_create( exodus_file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
275 if( NC_NOERR != fail )
276 {
277 reset_block( block_info );
278 return MB_FAILURE;
279 }
280
281 if( write_header( mesh_info, block_info, sideset_info, nodeset_info, exodus_file_name ) != MB_SUCCESS )
282 {
283 reset_block( block_info );
284 return MB_FAILURE;
285 }
286
287 {
288 // write dummy time_whole
289 double timev = 0.0; // dummy, to make paraview happy
290 size_t start = 0, count = 1;
291 int nc_var;
292 std::vector< int > dims;
293 GET_VAR( "time_whole", nc_var, dims );
294 fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &timev );
295 if( NC_NOERR != fail )
296 {
297 MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" );
298 }
299 }
300
301 if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
302 {
303 reset_block( block_info );
304 return MB_FAILURE;
305 }
306
307 if( !mesh_info.polyhedronFaces.empty() )
308 {
309 if( write_poly_faces( mesh_info ) != MB_SUCCESS )
310 {
311 reset_block( block_info );
312 return MB_FAILURE;
313 }
314 }
315
316 if( write_elementblocks( mesh_info, block_info ) )
317 {
318 reset_block( block_info );
319 return MB_FAILURE;
320 }
321
322 // Write the three maps
323 if( write_global_node_order_map( mesh_info.num_nodes, mesh_info.nodes ) != MB_SUCCESS )
324 {
325 reset_block( block_info );
326 return MB_FAILURE;
327 }
328
329 if( write_global_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
330 {
331 reset_block( block_info );
332 return MB_FAILURE;
333 }
334
335 if( write_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
336 {
337 reset_block( block_info );
338 return MB_FAILURE;
339 }
340
341 /*
342 if (write_elementmap(mesh_info) != MB_SUCCESS)
343 return MB_FAILURE;
344 */
345
346 if( write_BCs( sideset_info, nodeset_info ) != MB_SUCCESS )
347 {
348 reset_block( block_info );
349 return MB_FAILURE;
350 }
351
352 if( write_qa_records( mesh_info.qaRecords ) != MB_SUCCESS ) return MB_FAILURE;
353
354 // Copy the qa records into the argument
355 // mesh_info.qaRecords.swap(qa_records);
356 // Close the file
357 fail = nc_close( ncFile );
358 if( NC_NOERR != fail )
359 {
360 MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
361 }
362
363 return MB_SUCCESS;
364 }
365
366 ErrorCode WriteNCDF::gather_mesh_information( ExodusMeshInfo& mesh_info,
367 std::vector< MaterialSetData >& block_info,
368 std::vector< NeumannSetData >& sideset_info,
369 std::vector< DirichletSetData >& nodeset_info,
370 std::vector< EntityHandle >& blocks,
371 std::vector< EntityHandle >& sidesets,
372 std::vector< EntityHandle >& nodesets )
373 {
374 ErrorCode rval;
375 std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
376
377 mesh_info.num_nodes = 0;
378 mesh_info.num_elements = 0;
379 mesh_info.num_elementblocks = 0;
380 mesh_info.num_polyhedra_blocks = 0;
381
382 int id = 0;
383
384 vector_iter = blocks.begin();
385 end_vector_iter = blocks.end();
386
387 std::vector< EntityHandle > parent_meshsets;
388
389 // Clean out the bits for the element mark
390 rval = mdbImpl->tag_delete( mEntityMark );
391 if( MB_SUCCESS != rval ) return rval;
392 rval = mdbImpl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
393 if( MB_SUCCESS != rval ) return rval;
394
395 int highest_dimension_of_element_blocks = 0;
396
397 for( vector_iter = blocks.begin(); vector_iter != blocks.end(); ++vector_iter )
398 {
399 MaterialSetData block_data;
400
401 // For the purpose of qa records, get the parents of these blocks
402 if( mdbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
403
404 // Get all Entity Handles in the mesh set
405 Range dummy_range;
406 rval = mdbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
407 if( MB_SUCCESS != rval ) return rval;
408
409 // Skip empty blocks
410 if( dummy_range.empty() ) continue;
411
412 // Get the block's id
413 if( mdbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
414 {
415 MB_SET_ERR( MB_FAILURE, "Couldn't get block id from a tag for an element block" );
416 }
417
418 block_data.id = id;
419 block_data.number_attributes = 0;
420
421 // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
422
423 // Find the dimension of the last entity in this range
424 int this_dim = CN::Dimension( TYPE_FROM_HANDLE( dummy_range.back() ) );
425 if( this_dim > 3 )
426 {
427 MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "Block " << id << " contains entity sets" );
428 }
429 block_data.elements = dummy_range.subset_by_dimension( this_dim );
430
431 // End of -- wait a minute, we are doing some filtering here that doesn't make sense at this
432 // level CJS
433
434 // Get the entity type for this block, verifying that it's the same for all elements
435 EntityType entity_type = TYPE_FROM_HANDLE( block_data.elements.front() );
436 if( !block_data.elements.all_of_type( entity_type ) )
437 {
438 MB_SET_ERR( MB_FAILURE, "Entities in block " << id << " not of common type" );
439 }
440
441 int dimension = -1;
442 if( entity_type == MBQUAD || entity_type == MBTRI )
443 dimension = 2; // Output shells by default
444 else if( entity_type == MBEDGE )
445 dimension = 1;
446 else
447 dimension = CN::Dimension( entity_type );
448
449 if( dimension > highest_dimension_of_element_blocks ) highest_dimension_of_element_blocks = dimension;
450
451 std::vector< EntityHandle > tmp_conn;
452 rval = mdbImpl->get_connectivity( &( block_data.elements.front() ), 1, tmp_conn );
453 if( MB_SUCCESS != rval ) return rval;
454 block_data.element_type = ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
455
456 if( block_data.element_type == EXOII_MAX_ELEM_TYPE )
457 {
458 MB_SET_ERR( MB_FAILURE, "Element type in block " << id << " didn't get set correctly" );
459 }
460
461 if( block_data.element_type == EXOII_POLYGON )
462 {
463 // get all poly connectivity
464 int numconn = 0;
465 for( Range::iterator eit = block_data.elements.begin(); eit != block_data.elements.end(); eit++ )
466 {
467 EntityHandle polg = *eit;
468 int nnodes = 0;
469 const EntityHandle* conn = NULL;
470 rval = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
471 numconn += nnodes;
472 }
473 block_data.number_nodes_per_element = numconn;
474 }
475 else
476 block_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[block_data.element_type];
477
478 // Number of nodes for this block
479 block_data.number_elements = block_data.elements.size();
480
481 // Total number of elements
482 mesh_info.num_elements += block_data.number_elements;
483
484 // Get the nodes for the elements
485 rval = mWriteIface->gather_nodes_from_elements( block_data.elements, mEntityMark, mesh_info.nodes );
486 if( MB_SUCCESS != rval ) return rval;
487
488 // if polyhedra block
489 if( EXOII_POLYHEDRON == block_data.element_type )
490 {
491 rval = mdbImpl->get_connectivity( block_data.elements, mesh_info.polyhedronFaces );MB_CHK_ERR( rval );
492 mesh_info.num_polyhedra_blocks++;
493 }
494
495 if( !sidesets.empty() )
496 {
497 // If there are sidesets, keep track of which elements are being written out
498 for( Range::iterator iter = block_data.elements.begin(); iter != block_data.elements.end(); ++iter )
499 {
500 unsigned char bit = 0x1;
501 rval = mdbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
502 if( MB_SUCCESS != rval ) return rval;
503 }
504 }
505
506 block_info.push_back( block_data );
507
508 const void* data = NULL;
509 int size = 0;
510 if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mQaRecordTag, &( *vector_iter ), 1, &data, &size ) && NULL != data )
511 {
512 // There are qa records on this block - copy over to mesh qa records
513 const char* qa_rec = static_cast< const char* >( data );
514 int start = 0;
515 int count = 0;
516 for( int i = 0; i < size; i++ )
517 {
518 if( qa_rec[i] == '\0' )
519 {
520 std::string qa_string( &qa_rec[start], i - start );
521 mesh_info.qaRecords.push_back( qa_string );
522 start = i + 1;
523 count++;
524 }
525 }
526
527 // Constrained to multiples of 4 qa records
528 if( count > 0 ) assert( count % 4 == 0 );
529 }
530 }
531
532 mesh_info.num_elementblocks = block_info.size();
533
534 for( std::vector< MaterialSetData >::iterator blit = block_info.begin(); blit != block_info.end(); blit++ )
535 {
536 MaterialSetData& block = *blit;
537 if( block.element_type != EXOII_POLYHEDRON )
538 {
539 mesh_info.polyhedronFaces = subtract( mesh_info.polyhedronFaces, block.elements );
540 }
541 }
542
543 // If user hasn't entered dimension, we figure it out
544 if( mesh_info.num_dim == 0 )
545 {
546 // Never want 1 or zero dimensions
547 if( highest_dimension_of_element_blocks < 2 )
548 mesh_info.num_dim = 3;
549 else
550 mesh_info.num_dim = highest_dimension_of_element_blocks;
551 }
552
553 Range::iterator range_iter, end_range_iter;
554 range_iter = mesh_info.nodes.begin();
555 end_range_iter = mesh_info.nodes.end();
556
557 mesh_info.num_nodes = mesh_info.nodes.size();
558
559 //------nodesets--------
560
561 vector_iter = nodesets.begin();
562 end_vector_iter = nodesets.end();
563
564 for( ; vector_iter != end_vector_iter; ++vector_iter )
565 {
566 DirichletSetData nodeset_data;
567 nodeset_data.id = 0;
568 nodeset_data.number_nodes = 0;
569
570 // Get the nodeset's id
571 if( mdbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
572 {
573 MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for nodeset " << id );
574 }
575
576 nodeset_data.id = id;
577
578 std::vector< EntityHandle > node_vector;
579 // Get the nodes of the nodeset that are in mesh_info.nodes
580 if( mdbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
581 {
582 MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in nodeset " << id );
583 }
584
585 // Get the tag for distribution factors
586 const double* dist_factor_vector;
587 int dist_factor_size;
588 const void* ptr = 0;
589
590 int has_dist_factors = 0;
591 if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( *vector_iter ), 1, &ptr, &dist_factor_size ) == MB_SUCCESS &&
592 dist_factor_size )
593 has_dist_factors = 1;
594 dist_factor_size /= sizeof( double );
595 dist_factor_vector = reinterpret_cast< const double* >( ptr );
596 std::vector< EntityHandle >::iterator iter, end_iter;
597 iter = node_vector.begin();
598 end_iter = node_vector.end();
599
600 int j = 0;
601 unsigned char node_marked = 0;
602 ErrorCode result;
603 for( ; iter != end_iter; ++iter )
604 {
605 if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
606 result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
607
608 if( 0x1 == node_marked )
609 {
610 nodeset_data.nodes.push_back( *iter );
611 if( 0 != has_dist_factors )
612 nodeset_data.node_dist_factors.push_back( dist_factor_vector[j] );
613 else
614 nodeset_data.node_dist_factors.push_back( 1.0 );
615 }
616 j++;
617 }
618
619 nodeset_data.number_nodes = nodeset_data.nodes.size();
620 nodeset_info.push_back( nodeset_data );
621 }
622
623 //------sidesets--------
624 vector_iter = sidesets.begin();
625 end_vector_iter = sidesets.end();
626
627 for( ; vector_iter != end_vector_iter; ++vector_iter )
628 {
629 NeumannSetData sideset_data;
630
631 // Get the sideset's id
632 if( mdbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
633
634 sideset_data.id = id;
635 sideset_data.mesh_set_handle = *vector_iter;
636
637 // Get the sides in two lists, one forward the other reverse; starts with forward sense
638 // by convention
639 Range forward_elems, reverse_elems;
640 if( get_sideset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
641
642 ErrorCode result = get_valid_sides( forward_elems, mesh_info, 1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
643 result = get_valid_sides( reverse_elems, mesh_info, -1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
644
645 sideset_data.number_elements = sideset_data.elements.size();
646 sideset_info.push_back( sideset_data );
647 }
648
649 return MB_SUCCESS;
650 }
651
652 ErrorCode WriteNCDF::get_valid_sides( Range& elems,
653 ExodusMeshInfo& /*mesh_info*/,
654 const int sense,
655 NeumannSetData& sideset_data )
656 {
657 // This is where we see if underlying element of side set element is included in output
658
659 // Get the sideset-based info for distribution factors
660 const double* dist_factor_vector = 0;
661 int dist_factor_size = 0;
662
663 // Initialize dist_fac_iter to get rid of compiler warning
664 const double* dist_fac_iter = 0;
665 const void* ptr = 0;
666 bool has_dist_factors = false;
667 if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( sideset_data.mesh_set_handle ), 1, &ptr, &dist_factor_size ) ==
668 MB_SUCCESS &&
669 dist_factor_size )
670 {
671 has_dist_factors = true;
672 dist_factor_vector = reinterpret_cast< const double* >( ptr );
673 dist_fac_iter = dist_factor_vector;
674 dist_factor_size /= sizeof( double );
675 }
676
677 unsigned char element_marked = 0;
678 ErrorCode result;
679 for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
680 {
681 // Should insert here if "side" is a quad/tri on a quad/tri mesh
682 result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
683
684 if( 0x1 == element_marked )
685 {
686 sideset_data.elements.push_back( *iter );
687
688 // TJT TODO: the sense should really be # edges + 1or2
689 sideset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
690 }
691 else
692 { // Then "side" is probably a quad/tri on a hex/tet mesh
693 std::vector< EntityHandle > parents;
694 int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
695
696 // Get the adjacent parent element of "side"
697 if( mdbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
698 {
699 #if 0
700 // This is not treated as an error, print warning messages for
701 // debugging only
702 fprintf(stderr, "[Warning]: Couldn't get adjacencies for sideset.\n");
703 #endif
704 }
705
706 if( !parents.empty() )
707 {
708 // Make sure the adjacent parent element will be output
709 for( unsigned int k = 0; k < parents.size(); k++ )
710 {
711 result = mdbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
712
713 int side_no, this_sense, this_offset;
714 if( 0x1 == element_marked &&
715 mdbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
716 this_sense == sense )
717 {
718 sideset_data.elements.push_back( parents[k] );
719 sideset_data.side_numbers.push_back( side_no + 1 );
720 break;
721 }
722 }
723 }
724 else
725 {
726 #if 0
727 // This is not treated as an error, print warning messages for
728 // debugging only
729 fprintf(stderr, "[Warning]: No parent element exists for element in sideset %i\n", sideset_data.id);
730 #endif
731 }
732 }
733
734 if( sideset_data.elements.size() != 0 )
735 {
736 // Distribution factors
737 int num_nodes = CN::VerticesPerEntity( TYPE_FROM_HANDLE( *iter ) );
738 // put some dummy dist factors for polygons; why do we need them?
739 if( TYPE_FROM_HANDLE( *iter ) == MBPOLYGON ) num_nodes = 1; // dummy
740 if( has_dist_factors )
741 {
742 std::copy( dist_fac_iter, dist_fac_iter + num_nodes,
743 std::back_inserter( sideset_data.ss_dist_factors ) );
744 dist_fac_iter += num_nodes;
745 }
746 else
747 {
748 for( int j = 0; j < num_nodes; j++ )
749 sideset_data.ss_dist_factors.push_back( 1.0 );
750 }
751 }
752 }
753
754 return MB_SUCCESS;
755 }
756
757 ErrorCode WriteNCDF::write_qa_records( std::vector< std::string >& qa_record_list )
758 {
759 int i = 0;
760
761 for( std::vector< std::string >::iterator string_it = qa_record_list.begin(); string_it != qa_record_list.end(); )
762 {
763 for( int j = 0; j < 4; j++ )
764 write_qa_string( ( *string_it++ ).c_str(), i, j );
765 i++;
766 }
767
768 return MB_SUCCESS;
769 }
770
771 ErrorCode WriteNCDF::write_qa_string( const char* string, int record_number, int record_position )
772 {
773 // Get the variable id in the exodus file
774
775 std::vector< int > dims;
776 int temp_var = -1;
777 GET_VAR( "qa_records", temp_var, dims );
778 if( -1 == temp_var )
779 {
780 MB_SET_ERR( MB_FAILURE, "WriteNCDF:: Problem getting qa record variable" );
781 }
782 size_t count[3], start[3];
783
784 // Write out the record
785 start[0] = record_number;
786 start[1] = record_position;
787 start[2] = 0;
788
789 count[0] = 1;
790 count[1] = 1;
791 count[2] = (long)strlen( string ) + 1;
792 int fail = nc_put_vara_text( ncFile, temp_var, start, count, string );
793 if( NC_NOERR != fail )
794 {
795 MB_SET_ERR( MB_FAILURE, "Failed to position qa string variable" );
796 }
797
798 return MB_SUCCESS;
799 }
800
801 ErrorCode WriteNCDF::write_nodes( int num_nodes, Range& nodes, int dimension )
802 {
803 // Write coordinates names
804 int nc_var = -1;
805 std::vector< int > dims;
806 GET_VAR( "coor_names", nc_var, dims );
807 if( -1 == nc_var )
808 {
809 MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate name variable" );
810 }
811
812 size_t start[2] = { 0, 0 }, count[2] = { 1, ExoIIInterface::MAX_STR_LENGTH };
813 char dum_str[ExoIIInterface::MAX_STR_LENGTH];
814 strcpy( dum_str, "x" );
815 int fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
816 if( NC_NOERR != fail )
817 {
818 MB_SET_ERR( MB_FAILURE, "Trouble adding x coordinate name; netcdf message: " << nc_strerror( fail ) );
819 }
820
821 start[0] = 1;
822 strcpy( dum_str, "y" );
823 fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
824 if( NC_NOERR != fail )
825 {
826 MB_SET_ERR( MB_FAILURE, "Trouble adding y coordinate name; netcdf message: " << nc_strerror( fail ) );
827 }
828
829 start[0] = 2;
830 strcpy( dum_str, "z" );
831 fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
832 if( NC_NOERR != fail )
833 {
834 MB_SET_ERR( MB_FAILURE, "Trouble adding z coordinate name; netcdf message: " << nc_strerror( fail ) );
835 }
836
837 // See if should transform coordinates
838 ErrorCode result;
839 Tag trans_tag;
840 result = mdbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
841 bool transform_needed = true;
842 if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
843
844 int num_coords_to_fill = transform_needed ? 3 : dimension;
845
846 std::vector< double* > coord_arrays( 3 );
847 coord_arrays[0] = new double[num_nodes];
848 coord_arrays[1] = new double[num_nodes];
849 coord_arrays[2] = NULL;
850
851 if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
852
853 result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 1, coord_arrays );
854 if( result != MB_SUCCESS )
855 {
856 delete[] coord_arrays[0];
857 delete[] coord_arrays[1];
858 if( coord_arrays[2] ) delete[] coord_arrays[2];
859 return result;
860 }
861
862 if( transform_needed )
863 {
864 double trans_matrix[16];
865 const EntityHandle mesh = 0;
866 result = mdbImpl->tag_get_data( trans_tag, &mesh, 0, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
867
868 for( int i = 0; i < num_nodes; i++ )
869 {
870 double vec1[3];
871 double vec2[3];
872
873 vec2[0] = coord_arrays[0][i];
874 vec2[1] = coord_arrays[1][i];
875 vec2[2] = coord_arrays[2][i];
876
877 for( int row = 0; row < 3; row++ )
878 {
879 vec1[row] = 0.0;
880 for( int col = 0; col < 3; col++ )
881 vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
882 }
883
884 coord_arrays[0][i] = vec1[0];
885 coord_arrays[1][i] = vec1[1];
886 coord_arrays[2][i] = vec1[2];
887 }
888 }
889
890 // Write the nodes
891 nc_var = -1;
892 GET_VAR( "coord", nc_var, dims );
893 if( -1 == nc_var )
894 {
895 MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate variable" );
896 }
897 start[0] = 0;
898 count[1] = num_nodes;
899 fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[0][0] ) );
900 if( NC_NOERR != fail )
901 {
902 MB_SET_ERR( MB_FAILURE, "Trouble writing x coordinate" );
903 }
904
905 start[0] = 1;
906 fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[1][0] ) );
907 if( NC_NOERR != fail )
908 {
909 MB_SET_ERR( MB_FAILURE, "Trouble writing y coordinate" );
910 }
911
912 start[0] = 2;
913 fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[2][0] ) );
914 if( NC_NOERR != fail )
915 {
916 MB_SET_ERR( MB_FAILURE, "Trouble writing z coordinate" );
917 }
918
919 delete[] coord_arrays[0];
920 delete[] coord_arrays[1];
921 if( coord_arrays[2] ) delete[] coord_arrays[2];
922
923 return MB_SUCCESS;
924 }
925
926 ErrorCode WriteNCDF::write_poly_faces( ExodusMeshInfo& mesh_info )
927 {
928 // write all polygons that are not in another element block;
929 // usually they are nowhere else, but be sure, write in this block only ones that are not in the
930 // other blocks
931 Range pfaces = mesh_info.polyhedronFaces;
932
933 /*
934 * int fbconn1(num_nod_per_fa1) ;
935 fbconn1:elem_type = "nsided" ;
936 int fbepecnt1(num_fa_in_blk1) ;
937 fbepecnt1:entity_type1 = "NODE" ;
938 fbepecnt1:entity_type2 = "FACE" ;
939 */
940 if( pfaces.empty() ) return MB_SUCCESS;
941 char wname[CHAR_STR_LEN];
942 int nc_var = -1;
943 std::vector< int > dims;
944
945 // write one for each element block, to make paraview and visit happy
946 int num_faces_in_block = (int)pfaces.size();
947 for( unsigned int bl = 0; bl < mesh_info.num_polyhedra_blocks; bl++ )
948 {
949 INS_ID( wname, "fbconn%u", bl + 1 ); // it is the first block
950 GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
951
952 INS_ID( wname, "num_nod_per_fa%u", bl + 1 );
953 int ncdim, num_nod_per_face;
954 GET_DIM( ncdim, wname, num_nod_per_face );
955 int* connectivity = new int[num_nod_per_face];
956 int ixcon = 0, j = 0;
957 std::vector< int > fbepe( num_faces_in_block ); // fbepecnt1
958 for( Range::iterator eit = pfaces.begin(); eit != pfaces.end(); eit++ )
959 {
960 EntityHandle polyg = *eit;
961 int nnodes = 0;
962 const EntityHandle* conn = NULL;
963 ErrorCode rval = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
964 for( int k = 0; k < nnodes; k++ )
965 connectivity[ixcon++] = conn[k];
966 fbepe[j++] = nnodes;
967 }
968 size_t start[1] = { 0 }, count[1] = { 0 };
969 count[0] = ixcon;
970 int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
971 if( NC_NOERR != fail )
972 {
973 delete[] connectivity;
974 MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
975 }
976
977 INS_ID( wname, "fbepecnt%u", bl + 1 );
978 GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
979 count[0] = num_faces_in_block;
980
981 fail = nc_put_vara_int( ncFile, nc_var, start, count, &fbepe[0] );
982 if( NC_NOERR != fail )
983 {
984 delete[] connectivity;
985 MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
986 }
987
988 int id = bl + 1;
989 if( write_exodus_integer_variable( "fa_prop1", &id, bl, 1 ) != MB_SUCCESS )
990 {
991 MB_SET_ERR_CONT( "Problem writing element block id " << id );
992 }
993
994 int status = 1;
995 if( write_exodus_integer_variable( "fa_status", &status, bl, 1 ) != MB_SUCCESS )
996 {
997 MB_SET_ERR( MB_FAILURE, "Problem writing face block status" );
998 }
999
1000 delete[] connectivity;
1001 if( 0 == repeat_face_blocks ) break; // do not repeat face blocks
1002 }
1003
1004 return MB_SUCCESS;
1005 }
1006 ErrorCode WriteNCDF::write_header( ExodusMeshInfo& mesh_info,
1007 std::vector< MaterialSetData >& block_info,
1008 std::vector< NeumannSetData >& sideset_info,
1009 std::vector< DirichletSetData >& nodeset_info,
1010 const char* filename )
1011 {
1012 // Get the date and time
1013 char time[TIME_STR_LEN];
1014 char date[TIME_STR_LEN];
1015 time_and_date( time, date );
1016
1017 std::string title_string = "MOAB";
1018 title_string.append( "(" );
1019 title_string.append( filename );
1020 title_string.append( "): " );
1021 title_string.append( date );
1022 title_string.append( ": " );
1023 title_string.append( "time " );
1024
1025 if( title_string.length() > ExoIIInterface::MAX_LINE_LENGTH )
1026 title_string.resize( ExoIIInterface::MAX_LINE_LENGTH );
1027
1028 // Initialize the exodus file
1029
1030 int result = initialize_exodus_file( mesh_info, block_info, sideset_info, nodeset_info, title_string.c_str() );
1031
1032 if( result == MB_FAILURE ) return MB_FAILURE;
1033
1034 return MB_SUCCESS;
1035 }
1036
1037 ErrorCode WriteNCDF::write_elementblocks( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data )
1038 {
1039 unsigned int i;
1040 int block_index = 0; // Index into block list, may != 1 if there are inactive blocks
1041 int exodus_id = 1;
1042
1043 for( i = 0; i < block_data.size(); i++ )
1044 {
1045 MaterialSetData& block = block_data[i];
1046
1047 unsigned int num_nodes_per_elem = block.number_nodes_per_element;
1048
1049 // Write out the id for the block
1050
1051 int id = block.id;
1052 int num_values = 1;
1053
1054 if( write_exodus_integer_variable( "eb_prop1", &id, block_index, num_values ) != MB_SUCCESS )
1055 {
1056 MB_SET_ERR_CONT( "Problem writing element block id " << id );
1057 }
1058
1059 // Write out the block status
1060
1061 int status = 1;
1062 if( 0 == block.number_elements )
1063 {
1064 MB_SET_ERR( MB_FAILURE, "No elements in block " << id );
1065 }
1066
1067 if( write_exodus_integer_variable( "eb_status", &status, block_index, num_values ) != MB_SUCCESS )
1068 {
1069 MB_SET_ERR( MB_FAILURE, "Problem writing element block status" );
1070 }
1071
1072 //
1073 // Map the connectivity to the new nodes
1074 const unsigned int num_elem = block.number_elements;
1075 unsigned int num_nodes = num_nodes_per_elem * num_elem;
1076 if( EXOII_POLYGON == block.element_type || EXOII_POLYHEDRON == block.element_type )
1077 {
1078 num_nodes = num_nodes_per_elem;
1079 }
1080 int* connectivity = new int[num_nodes];
1081
1082 ErrorCode result = MB_SUCCESS;
1083 if( block.element_type != EXOII_POLYHEDRON )
1084 mWriteIface->get_element_connect( num_elem, num_nodes_per_elem, mGlobalIdTag, block.elements, mGlobalIdTag,
1085 exodus_id, connectivity );
1086 if( result != MB_SUCCESS )
1087 {
1088 delete[] connectivity;
1089 MB_SET_ERR( result, "Couldn't get element array to write from" );
1090 }
1091
1092 // If necessary, convert from EXODUS to CN node order
1093 const EntityType elem_type = ExoIIUtil::ExoIIElementMBEntity[block.element_type];
1094 assert( block.elements.all_of_type( elem_type ) );
1095 const int* reorder = 0;
1096 if( block.element_type != EXOII_POLYHEDRON && block.element_type != EXOII_POLYGON )
1097 reorder = exodus_elem_order_map[elem_type][block.number_nodes_per_element];
1098 if( reorder )
1099 WriteUtilIface::reorder( reorder, connectivity, block.number_elements, block.number_nodes_per_element );
1100
1101 char wname[CHAR_STR_LEN];
1102 int nc_var = -1;
1103 std::vector< int > dims;
1104 if( block.element_type != EXOII_POLYHEDRON )
1105 {
1106 exodus_id += num_elem;
1107 INS_ID( wname, "connect%u", i + 1 );
1108
1109 GET_VAR( wname, nc_var, dims );
1110 if( -1 == nc_var )
1111 {
1112 delete[] connectivity;
1113 MB_SET_ERR( MB_FAILURE, "Couldn't get connectivity variable" );
1114 }
1115 }
1116
1117 if( EXOII_POLYGON == block.element_type )
1118 {
1119 size_t start[1] = { 0 }, count[1] = { num_nodes_per_elem };
1120 int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1121 if( NC_NOERR != fail )
1122 {
1123 delete[] connectivity;
1124 MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
1125 }
1126 // now put also number ebepecnt1
1127 INS_ID( wname, "ebepecnt%u", i + 1 );
1128 GET_VAR( wname, nc_var, dims );
1129 count[0] = block.number_elements;
1130 start[0] = 0;
1131 // reuse connectivity array, to not allocate another one
1132 int j = 0;
1133 for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); j++, eit++ )
1134 {
1135 EntityHandle polg = *eit;
1136 int nnodes = 0;
1137 const EntityHandle* conn = NULL;
1138 ErrorCode rval = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
1139 connectivity[j] = nnodes;
1140 }
1141 fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1142 if( NC_NOERR != fail )
1143 {
1144 delete[] connectivity;
1145 MB_SET_ERR( MB_FAILURE, "Couldn't write ebepecnt variable" );
1146 }
1147 }
1148 else if( block.element_type != EXOII_POLYHEDRON )
1149 {
1150 size_t start[2] = { 0, 0 }, count[2] = { num_elem, num_nodes_per_elem };
1151 int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1152 if( NC_NOERR != fail )
1153 {
1154 delete[] connectivity;
1155 MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
1156 }
1157 }
1158 else // if (block.element_type == EXOII_POLYHEDRON)
1159 {
1160 /* write a lot of stuff // faconn
1161 num_fa_in_blk1 = 15 ;
1162 num_nod_per_fa1 = 58 ;
1163 num_el_in_blk1 = 3 ;
1164 num_fac_per_el1 = 17 ;
1165 int fbconn1(num_nod_per_fa1) ;
1166 fbconn1:elem_type = "nsided" ;
1167 int fbepecnt1(num_fa_in_blk1) ;
1168 fbepecnt1:entity_type1 = "NODE" ;
1169 fbepecnt1:entity_type2 = "FACE" ;
1170 int facconn1(num_fac_per_el1) ;
1171 facconn1:elem_type = "NFACED" ;
1172 int ebepecnt1(num_el_in_blk1) ;
1173 ebepecnt1:entity_type1 = "FACE" ;
1174 ebepecnt1:entity_type2 = "ELEM" ;
1175
1176 */
1177 Range& block_faces = mesh_info.polyhedronFaces;
1178 // ErrorCode rval = mdbImpl->get_connectivity(block.elements, block_faces);
1179 // MB_CHK_ERR(rval);
1180
1181 // reuse now connectivity for facconn1
1182 INS_ID( wname, "facconn%u", i + 1 );
1183 GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
1184
1185 std::vector< int > ebepe( block.elements.size() ); // ebepecnt1
1186 int ixcon = 0, j = 0;
1187 size_t start[1] = { 0 }, count[1] = { 0 };
1188
1189 for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
1190 {
1191 EntityHandle polyh = *eit;
1192 int nfaces = 0;
1193 const EntityHandle* conn = NULL;
1194 ErrorCode rval = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
1195 for( int k = 0; k < nfaces; k++ )
1196 {
1197 int index = block_faces.index( conn[k] );
1198 if( index == -1 ) MB_SET_ERR( MB_FAILURE, "Couldn't find face in polyhedron" );
1199 connectivity[ixcon++] = index + 1;
1200 }
1201 ebepe[j++] = nfaces;
1202 // num_faces+=nfaces;
1203 }
1204 count[0] = ixcon; // facconn1
1205 int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1206 if( NC_NOERR != fail )
1207 {
1208 delete[] connectivity;
1209 MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
1210 }
1211
1212 INS_ID( wname, "ebepecnt%u", i + 1 );
1213 GET_VAR( wname, nc_var, dims ); // ebepecnt# variable, 1 dimensional
1214 count[0] = block.elements.size();
1215
1216 fail = nc_put_vara_int( ncFile, nc_var, start, count, &ebepe[0] );
1217 if( NC_NOERR != fail )
1218 {
1219 delete[] connectivity;
1220 MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
1221 }
1222 }
1223 block_index++;
1224 delete[] connectivity;
1225 }
1226
1227 return MB_SUCCESS;
1228 }
1229
1230 ErrorCode WriteNCDF::write_global_node_order_map( int num_nodes, Range& nodes )
1231 {
1232 // Note: this routine bypasses the standard exodusII interface for efficiency!
1233
1234 // Node order map
1235 int* map = new int[num_nodes];
1236
1237 // For now, output a dummy map!
1238
1239 Range::iterator range_iter, end_iter;
1240 range_iter = nodes.begin();
1241 end_iter = nodes.end();
1242
1243 int i = 0;
1244
1245 for( ; range_iter != end_iter; ++range_iter )
1246 {
1247 // TODO -- do we really want to cast this to an int?
1248 map[i++] = (int)ID_FROM_HANDLE( *range_iter );
1249 }
1250
1251 // Output array and cleanup
1252
1253 int error = write_exodus_integer_variable( "node_num_map", map, 0, num_nodes );
1254
1255 if( map ) delete[] map;
1256
1257 if( error < 0 )
1258 {
1259 MB_SET_ERR( MB_FAILURE, "Failed writing global node order map" );
1260 }
1261
1262 return MB_SUCCESS;
1263 }
1264
1265 ErrorCode WriteNCDF::write_global_element_order_map( int num_elements )
1266 {
1267 // Allocate map array
1268 int* map = new int[num_elements];
1269
1270 // Many Sandia codes assume this map is unique, and CUBIT does not currently
1271 // have unique ids for all elements. Therefore, to make sure nothing crashes,
1272 // insert a dummy map...
1273
1274 for( int i = 0; i < num_elements; i++ )
1275 map[i] = i + 1;
1276
1277 // Output array and cleanup
1278
1279 int error = write_exodus_integer_variable( "elem_num_map", map, 0, num_elements );
1280
1281 if( map ) delete[] map;
1282
1283 if( error < 0 )
1284 {
1285 MB_SET_ERR( MB_FAILURE, "Failed writing global element order map" );
1286 }
1287
1288 return MB_SUCCESS;
1289 }
1290
1291 ErrorCode WriteNCDF::write_element_order_map( int num_elements )
1292 {
1293 // Note: this routine bypasses the standard exodusII interface for efficiency!
1294
1295 // Element order map
1296 int* map = new int[num_elements];
1297
1298 // For now, output a dummy map!
1299
1300 for( int i = 0; i < num_elements; i++ )
1301 {
1302 map[i] = i + 1;
1303 }
1304
1305 // Output array and cleanup
1306
1307 int error = write_exodus_integer_variable( "elem_map", map, 0, num_elements );
1308
1309 if( map ) delete[] map;
1310
1311 if( error < 0 )
1312 {
1313 MB_SET_ERR( MB_FAILURE, "Failed writing element map" );
1314 }
1315
1316 return MB_SUCCESS;
1317 }
1318
1319 ErrorCode WriteNCDF::write_exodus_integer_variable( const char* variable_name,
1320 int* variable_array,
1321 int start_position,
1322 int number_values )
1323 {
1324 // Note: this routine bypasses the standard exodusII interface for efficiency!
1325
1326 // Write directly to netcdf interface for efficiency
1327
1328 // Get the variable id of the element map
1329 int nc_var = -1;
1330 std::vector< int > dims;
1331 GET_VAR( variable_name, nc_var, dims );
1332 if( -1 == nc_var )
1333 {
1334 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate variable " << variable_name << " in file" );
1335 }
1336 // This contortion is necessary because netCDF is expecting nclongs;
1337 // fortunately it's necessary only when ints and nclongs aren't the same size
1338
1339 size_t start[1], count[1];
1340 start[0] = start_position;
1341 count[0] = number_values;
1342
1343 int fail = NC_NOERR;
1344 if( sizeof( int ) == sizeof( long ) )
1345 {
1346 fail = nc_put_vara_int( ncFile, nc_var, start, count, variable_array );
1347 }
1348 else
1349 {
1350 long* lptr = new long[number_values];
1351 for( int jj = 0; jj < number_values; jj++ )
1352 lptr[jj] = variable_array[jj];
1353 fail = nc_put_vara_long( ncFile, nc_var, start, count, lptr );
1354 delete[] lptr;
1355 }
1356
1357 if( NC_NOERR != fail )
1358 {
1359 MB_SET_ERR( MB_FAILURE, "Failed to store variable " << variable_name );
1360 }
1361
1362 return MB_SUCCESS;
1363 }
1364
1365 ErrorCode WriteNCDF::write_BCs( std::vector< NeumannSetData >& sidesets, std::vector< DirichletSetData >& nodesets )
1366 {
1367 unsigned int i, j;
1368 int id;
1369 int ns_index = -1;
1370
1371 for( std::vector< DirichletSetData >::iterator ns_it = nodesets.begin(); ns_it != nodesets.end(); ++ns_it )
1372 {
1373 // Get number of nodes in set
1374 int number_nodes = ( *ns_it ).number_nodes;
1375 if( 0 == number_nodes ) continue;
1376
1377 // If we're here, we have a non-empty nodeset; increment the index
1378 ns_index++;
1379
1380 // Get the node set id
1381 id = ( *ns_it ).id;
1382
1383 // Build new array to old exodus ids
1384 int* exodus_id_array = new int[number_nodes];
1385 double* dist_factor_array = new double[number_nodes];
1386
1387 std::vector< EntityHandle >::iterator begin_iter, end_iter;
1388 std::vector< double >::iterator other_iter;
1389 begin_iter = ( *ns_it ).nodes.begin();
1390 end_iter = ( *ns_it ).nodes.end();
1391 other_iter = ( *ns_it ).node_dist_factors.begin();
1392
1393 j = 0;
1394 int exodus_id;
1395 ErrorCode result;
1396 // Fill up node array and dist. factor array at the same time
1397 for( ; begin_iter != end_iter; ++begin_iter )
1398 {
1399 result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting id tag data" );
1400
1401 exodus_id_array[j] = exodus_id;
1402 dist_factor_array[j] = *( other_iter );
1403 ++other_iter;
1404 j++;
1405 }
1406
1407 // Write out the id for the nodeset
1408
1409 int num_values = 1;
1410
1411 result = write_exodus_integer_variable( "ns_prop1", &id, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
1412
1413 // Write out the nodeset status
1414
1415 int status = 1;
1416 if( !number_nodes ) status = 0;
1417
1418 result = write_exodus_integer_variable( "ns_status", &status, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set status", MB_FAILURE );
1419
1420 // Write it out
1421 char wname[CHAR_STR_LEN];
1422 int nc_var = -1;
1423 std::vector< int > dims;
1424 INS_ID( wname, "node_ns%d", ns_index + 1 );
1425 GET_VAR( wname, nc_var, dims );
1426 if( -1 == nc_var )
1427 {
1428 MB_SET_ERR( MB_FAILURE, "Failed to get node_ns variable" );
1429 }
1430
1431 size_t start = 0, count = number_nodes;
1432 int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, exodus_id_array );
1433 if( NC_NOERR != fail )
1434 {
1435 MB_SET_ERR( MB_FAILURE, "Failed writing exodus id array" );
1436 }
1437
1438 // Write out nodeset distribution factors
1439 INS_ID( wname, "dist_fact_ns%d", ns_index + 1 );
1440 nc_var = -1;
1441 GET_VAR( wname, nc_var, dims );
1442 if( -1 == nc_var )
1443 {
1444 MB_SET_ERR( MB_FAILURE, "Failed to get dist_fact variable" );
1445 }
1446 fail = nc_put_vara_double( ncFile, nc_var, &start, &count, dist_factor_array );
1447 if( NC_NOERR != fail )
1448 {
1449 MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" );
1450 }
1451
1452 delete[] dist_factor_array;
1453 delete[] exodus_id_array;
1454 }
1455
1456 // Now do sidesets
1457 int ss_index = 0; // Index of sideset - not the same as 'i' because
1458 // only writing non-empty side sets
1459 for( i = 0; i < sidesets.size(); i++ )
1460 {
1461 NeumannSetData sideset_data = sidesets[i];
1462
1463 // Get the side set id
1464 int side_set_id = sideset_data.id;
1465
1466 // Get number of elements in set
1467 int number_elements = sideset_data.number_elements;
1468 if( 0 == number_elements ) continue;
1469
1470 // Build new array to old exodus ids
1471 int* output_element_ids = new int[number_elements];
1472 int* output_element_side_numbers = new int[number_elements];
1473
1474 std::vector< EntityHandle >::iterator begin_iter, end_iter;
1475 begin_iter = sideset_data.elements.begin();
1476 end_iter = sideset_data.elements.end();
1477 std::vector< int >::iterator side_iter = sideset_data.side_numbers.begin();
1478
1479 // Get the tag handle
1480 j = 0;
1481 int exodus_id;
1482
1483 // For each "side"
1484 for( ; begin_iter != end_iter; ++begin_iter, ++side_iter )
1485 {
1486 ErrorCode result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting exodus id for sideset element "
1487 << (long unsigned int)ID_FROM_HANDLE( *begin_iter ) );
1488
1489 output_element_ids[j] = exodus_id;
1490 output_element_side_numbers[j++] = *side_iter;
1491 }
1492
1493 if( 0 != number_elements )
1494 {
1495 // Write out the id for the nodeset
1496
1497 int num_values = 1;
1498
1499 // ss_prop1[ss_index] = side_set_id
1500 ErrorCode result = write_exodus_integer_variable( "ss_prop1", &side_set_id, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
1501
1502 // FIXME : Something seems wrong here. The we are within a block
1503 // started with if (0 != number_elements), so this condition is always
1504 // false. This code seems to indicate that we want to write all
1505 // sidesets, not just empty ones. But the code both here and in
1506 // initialize_exodus_file() skip empty side sets.
1507 // - j.k. 2007-03-09
1508 int status = 1;
1509 if( 0 == number_elements ) status = 0;
1510
1511 // ss_status[ss_index] = status
1512 result = write_exodus_integer_variable( "ss_status", &status, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing side set status", MB_FAILURE );
1513
1514 // Increment ss_index now because we want a) we need to
1515 // increment it somewhere within the if (0 != number_elements)
1516 // block and b) the above calls need a zero-based index
1517 // while the following use a one-based index.
1518 ++ss_index;
1519
1520 char wname[CHAR_STR_LEN];
1521 int nc_var;
1522 std::vector< int > dims;
1523 INS_ID( wname, "elem_ss%d", ss_index );
1524 GET_VAR( wname, nc_var, dims );
1525 if( -1 == nc_var )
1526 {
1527 MB_SET_ERR( MB_FAILURE, "Failed to get elem_ss variable" );
1528 }
1529 size_t start = 0, count = number_elements;
1530 int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_ids );
1531 if( NC_NOERR != fail )
1532 {
1533 MB_SET_ERR( MB_FAILURE, "Failed writing sideset element array" );
1534 }
1535
1536 INS_ID( wname, "side_ss%d", ss_index );
1537 nc_var = -1;
1538 GET_VAR( wname, nc_var, dims );
1539 if( -1 == nc_var )
1540 {
1541 MB_SET_ERR( MB_FAILURE, "Failed to get side_ss variable" );
1542 }
1543 fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_side_numbers );
1544 if( NC_NOERR != fail )
1545 {
1546 MB_SET_ERR( MB_FAILURE, "Failed writing sideset side array" );
1547 }
1548
1549 INS_ID( wname, "dist_fact_ss%d", ss_index );
1550 nc_var = -1;
1551 GET_VAR( wname, nc_var, dims );
1552 if( -1 == nc_var )
1553 {
1554 MB_SET_ERR( MB_FAILURE, "Failed to get sideset dist factors variable" );
1555 }
1556 count = sideset_data.ss_dist_factors.size();
1557 fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &( sideset_data.ss_dist_factors[0] ) );
1558 if( NC_NOERR != fail )
1559 {
1560 MB_SET_ERR( MB_FAILURE, "Failed writing sideset dist factors array" );
1561 }
1562 }
1563
1564 delete[] output_element_ids;
1565 delete[] output_element_side_numbers;
1566 }
1567
1568 return MB_SUCCESS;
1569 }
1570
1571 ErrorCode WriteNCDF::initialize_exodus_file( ExodusMeshInfo& mesh_info,
1572 std::vector< MaterialSetData >& block_data,
1573 std::vector< NeumannSetData >& sideset_data,
1574 std::vector< DirichletSetData >& nodeset_data,
1575 const char* title_string,
1576 bool write_maps,
1577 bool /* write_sideset_distribution_factors */ )
1578 {
1579 // This routine takes the place of the exodusii routine ex_put_init,
1580 // and additionally pre-defines variables such as qa, element blocks,
1581 // sidesets and nodesets in a single pass.
1582 //
1583 // This is necessary because of the way exodusII works. Each time the
1584 // netcdf routine endef is called to take the file out of define mode,
1585 // the entire file is copied before the new information is added.
1586 //
1587 // With very large files, this is simply not workable. This routine takes
1588 // the definition portions of all applicable exodus routines and puts them
1589 // in a single definition, preventing repeated copying of the file.
1590 //
1591 // Most of the code is copied directly from the applicable exodus routine,
1592 // and thus this routine may not seem very consistent in usage or variable
1593 // naming, etc.
1594
1595 // Perform the initializations
1596
1597 int element_block_index;
1598
1599 // Inquire on defined string dimension and general dimension for qa
1600
1601 int dim_str, dim_four, dim_line, dim_time;
1602 if( nc_def_dim( ncFile, "len_string", ExoIIInterface::MAX_STR_LENGTH, &dim_str ) != NC_NOERR )
1603 {
1604 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get string length in file" );
1605 }
1606
1607 if( nc_def_dim( ncFile, "len_line", ExoIIInterface::MAX_STR_LENGTH, &dim_line ) != NC_NOERR )
1608 {
1609 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get line length in file" );
1610 }
1611
1612 if( nc_def_dim( ncFile, "four", 4, &dim_four ) != NC_NOERR )
1613 {
1614 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate four in file" );
1615 }
1616
1617 if( nc_def_dim( ncFile, "time_step", 1, &dim_time ) != NC_NOERR )
1618 {
1619 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate time step in file" );
1620 }
1621 // some whole_time dummy :(
1622 int dtime;
1623 if( NC_NOERR != nc_def_var( ncFile, "time_whole", NC_DOUBLE, 1, &dim_time, &dtime ) )
1624 {
1625 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define time whole array" );
1626 }
1627
1628 /* Put file into define mode */
1629
1630 // It is possible that an NT filename using backslashes is in the title string
1631 // this can give fits to unix codes where the backslash is assumed to be an escape
1632 // sequence. For the exodus file, just flip the backslash to a slash to prevent
1633 // this problem
1634
1635 // Get a working copy of the title_string;
1636
1637 char working_title[80];
1638 strncpy( working_title, title_string, 79 );
1639
1640 int length = strlen( working_title );
1641 for( int pos = 0; pos < length; pos++ )
1642 {
1643 if( working_title[pos] == '\\' ) working_title[pos] = '/';
1644 }
1645
1646 if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "title", length, working_title ) )
1647 {
1648 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define title attribute" );
1649 }
1650
1651 // Add other attributes while we're at it
1652 float dum_vers = 6.28F;
1653 if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "api_version", NC_FLOAT, 1, &dum_vers ) )
1654 {
1655 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define api_version attribute" );
1656 }
1657 dum_vers = 6.28F;
1658 if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "version", NC_FLOAT, 1, &dum_vers ) )
1659 {
1660 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define version attribute" );
1661 }
1662 int dum_siz = sizeof( double );
1663 if( NC_NOERR != nc_put_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &dum_siz ) )
1664 {
1665 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define floating pt word size attribute" );
1666 }
1667
1668 // Set up number of dimensions
1669
1670 int num_el_blk, num_elem, num_nodes, num_dim, num_fa_blk, num_faces;
1671 if( nc_def_dim( ncFile, "num_dim", (size_t)mesh_info.num_dim, &num_dim ) != NC_NOERR )
1672 {
1673 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dimensions" );
1674 }
1675
1676 if( nc_def_dim( ncFile, "num_nodes", mesh_info.num_nodes, &num_nodes ) != NC_NOERR )
1677 {
1678 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
1679 }
1680
1681 int num_nod_per_fa; // it is needed for polyhedron only; need to compute it (connectivity of
1682 // faces of polyhedra)
1683 if( mesh_info.polyhedronFaces.size() > 0 )
1684 if( nc_def_dim( ncFile, "num_faces", (int)mesh_info.polyhedronFaces.size(), &num_faces ) != NC_NOERR )
1685 {
1686 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
1687 }
1688
1689 if( nc_def_dim( ncFile, "num_elem", mesh_info.num_elements, &num_elem ) != NC_NOERR )
1690 {
1691 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements" );
1692 }
1693
1694 if( nc_def_dim( ncFile, "num_el_blk", mesh_info.num_elementblocks, &num_el_blk ) != NC_NOERR )
1695 {
1696 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of element blocks" );
1697 }
1698
1699 /* ...and some variables */
1700
1701 /* Element block id status array */
1702 int idstat = -1;
1703 if( NC_NOERR != nc_def_var( ncFile, "eb_status", NC_LONG, 1, &num_el_blk, &idstat ) )
1704 {
1705 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block status array" );
1706 }
1707
1708 /* Element block id array */
1709
1710 int idarr = -1;
1711 if( NC_NOERR != nc_def_var( ncFile, "eb_prop1", NC_LONG, 1, &num_el_blk, &idarr ) )
1712 {
1713 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block id array" );
1714 }
1715
1716 /* store property name as attribute of property array variable */
1717 if( NC_NOERR != nc_put_att_text( ncFile, idarr, "name", strlen( "ID" ), "ID" ) )
1718 {
1719 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element block property name ID" );
1720 }
1721
1722 // count how many are polyhedron blocks
1723 int num_fa_blocks = 0; //, num_polyh_blocks = 0;
1724 for( unsigned int i = 0; i < block_data.size(); i++ )
1725 {
1726 MaterialSetData& block = block_data[i];
1727 if( EXOII_POLYHEDRON == block.element_type )
1728 {
1729 num_fa_blocks++;
1730 // num_polyh_blocks++;
1731 }
1732 }
1733 if( 0 == this->repeat_face_blocks && num_fa_blocks > 1 ) num_fa_blocks = 1;
1734 char wname[CHAR_STR_LEN];
1735
1736 if( num_fa_blocks > 0 )
1737 {
1738 /* face block id status array */
1739 if( nc_def_dim( ncFile, "num_fa_blk", num_fa_blocks, &num_fa_blk ) != NC_NOERR )
1740 {
1741 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of face blocks" );
1742 }
1743
1744 int idstatf = -1;
1745 if( NC_NOERR != nc_def_var( ncFile, "fa_status", NC_LONG, 1, &num_fa_blk, &idstatf ) )
1746 {
1747 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block status array" );
1748 }
1749
1750 /* Element block id array */
1751
1752 int idarrf = -1;
1753 if( NC_NOERR != nc_def_var( ncFile, "fa_prop1", NC_LONG, 1, &num_fa_blk, &idarrf ) )
1754 {
1755 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block id array" );
1756 }
1757
1758 /* store property name as attribute of property array variable */
1759 if( NC_NOERR != nc_put_att_text( ncFile, idarrf, "name", strlen( "ID" ), "ID" ) )
1760 {
1761 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store face block property name ID" );
1762 }
1763 // determine the number of num_nod_per_face
1764 /*
1765 num_fa_in_blk1 = 15 ;
1766 num_nod_per_fa1 = 58 ;
1767
1768 int fbconn1(num_nod_per_fa1) ;
1769 fbconn1:elem_type = "nsided" ;
1770 int fbepecnt1(num_fa_in_blk1) ;
1771 fbepecnt1:entity_type1 = "NODE" ;
1772 fbepecnt1:entity_type2 = "FACE" ;
1773 */
1774
1775 int num_nodes_per_face = 0;
1776
1777 int dims[1]; // maybe 1 is enough here
1778 for( Range::iterator eit = mesh_info.polyhedronFaces.begin(); eit != mesh_info.polyhedronFaces.end(); eit++ )
1779 {
1780 EntityHandle polyg = *eit;
1781 int nnodes = 0;
1782 const EntityHandle* conn = NULL;
1783 ErrorCode rval = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
1784 num_nodes_per_face += nnodes;
1785 }
1786
1787 // duplicate if needed; default is not duplicate
1788 for( int j = 1; j <= num_fa_blocks; j++ )
1789 {
1790 INS_ID( wname, "num_nod_per_fa%d", j );
1791 if( nc_def_dim( ncFile, wname, (size_t)num_nodes_per_face, &num_nod_per_fa ) != NC_NOERR )
1792 {
1793 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
1794 }
1795 dims[0] = num_nod_per_fa;
1796 INS_ID( wname, "fbconn%d", j ); // first one, or more
1797 int fbconn;
1798 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbconn ) )
1799 {
1800 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for face block " << 1 );
1801 }
1802 std::string element_type_string( "nsided" );
1803 if( NC_NOERR != nc_put_att_text( ncFile, fbconn, "elem_type", element_type_string.length(),
1804 element_type_string.c_str() ) )
1805 {
1806 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type nsided " );
1807 }
1808
1809 INS_ID( wname, "num_fa_in_blk%d", j ); // first one, or more
1810 int num_fa_in_blk;
1811 if( nc_def_dim( ncFile, wname, (size_t)mesh_info.polyhedronFaces.size(), &num_fa_in_blk ) != NC_NOERR )
1812 {
1813 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
1814 }
1815
1816 // fbepecnt
1817 INS_ID( wname, "fbepecnt%d", j ); // first one, or more
1818 int fbepecnt;
1819 dims[0] = num_fa_in_blk;
1820 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbepecnt ) )
1821 {
1822 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create fbepecnt array for block " << 1 );
1823 }
1824 std::string enttype1( "NODE" );
1825 if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type1", enttype1.length(), enttype1.c_str() ) )
1826 {
1827 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 1 " );
1828 }
1829 std::string enttype2( "FACE" );
1830 if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type2", enttype2.length(), enttype2.c_str() ) )
1831 {
1832 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 2 " );
1833 }
1834 }
1835 }
1836
1837 // Define element blocks
1838
1839 for( unsigned int i = 0; i < block_data.size(); i++ )
1840 {
1841 MaterialSetData& block = block_data[i];
1842
1843 element_block_index = i + 1;
1844 int num_el_in_blk = -1, num_att_in_blk = -1;
1845 int blk_attrib, connect;
1846
1847 /* Define number of elements in this block */
1848
1849 INS_ID( wname, "num_el_in_blk%d", element_block_index );
1850 if( nc_def_dim( ncFile, wname, (size_t)block.number_elements, &num_el_in_blk ) != NC_NOERR )
1851 {
1852 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements/block for block " << i + 1 );
1853 }
1854
1855 /* Define number of nodes per element for this block */
1856 INS_ID( wname, "num_nod_per_el%d", element_block_index );
1857 int num_nod_per_el = -1;
1858 if( EXOII_POLYHEDRON != block.element_type )
1859 if( nc_def_dim( ncFile, wname, (size_t)block.number_nodes_per_element, &num_nod_per_el ) != NC_NOERR )
1860 {
1861 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes/element for block " << block.id );
1862 }
1863
1864 /* Define element attribute array for this block */
1865 int dims[3];
1866 if( block.number_attributes > 0 )
1867 {
1868 INS_ID( wname, "num_att_in_blk%d", element_block_index );
1869 if( nc_def_dim( ncFile, wname, (size_t)block.number_attributes, &num_att_in_blk ) != NC_NOERR )
1870 {
1871 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of attributes in block " << block.id );
1872 }
1873
1874 INS_ID( wname, "attrib%d", element_block_index );
1875 dims[0] = num_el_in_blk;
1876 dims[1] = num_att_in_blk;
1877 if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 2, dims, &blk_attrib ) )
1878 {
1879 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define attributes for element block " << block.id );
1880 }
1881 }
1882
1883 /* Define element connectivity array for this block */
1884
1885 if( EXOII_POLYGON != block.element_type && EXOII_POLYHEDRON != block.element_type )
1886 {
1887 INS_ID( wname, "connect%d", element_block_index );
1888 dims[0] = num_el_in_blk;
1889 dims[1] = num_nod_per_el;
1890 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 2, dims, &connect ) )
1891 {
1892 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
1893 }
1894
1895 /* Store element type as attribute of connectivity variable */
1896
1897 std::string element_type_string( ExoIIUtil::ElementTypeNames[block.element_type] );
1898 if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1899 element_type_string.c_str() ) )
1900 {
1901 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
1902 }
1903 }
1904 else if( EXOII_POLYGON == block.element_type )
1905 {
1906 INS_ID( wname, "connect%d", element_block_index );
1907 // need to define num_nod_per_el as total number of nodes
1908 // ebepecnt1 as number of nodes per polygon
1909 /*
1910 * int connect1(num_nod_per_el1) ;
1911 connect1:elem_type = "nsided" ;
1912 int ebepecnt1(num_el_in_blk1) ;
1913 ebepecnt1:entity_type1 = "NODE" ;
1914 ebepecnt1:entity_type2 = "ELEM" ;*/
1915 dims[0] = num_nod_per_el;
1916 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &connect ) )
1917 {
1918 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
1919 }
1920 std::string element_type_string( "nsided" );
1921 if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1922 element_type_string.c_str() ) )
1923 {
1924 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
1925 }
1926 INS_ID( wname, "ebepecnt%d", element_block_index );
1927 int ebepecnt;
1928 dims[0] = num_el_in_blk;
1929 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
1930 {
1931 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
1932 }
1933 std::string etype1( "NODE" );
1934 if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
1935 {
1936 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
1937 }
1938 std::string etype2( "ELEM" );
1939 if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
1940 {
1941 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
1942 }
1943 }
1944 else if( EXOII_POLYHEDRON == block.element_type )
1945 {
1946 // INS_ID(wname, "connect%d", element_block_index);
1947 /*
1948 testn face example: 3 polyh, 15 total faces, 2 shared; 15+2 = 17
1949 num_elem = 3 ;
1950 num_face = 15 ; // not needed?
1951 num_el_blk = 1 ;
1952
1953 num_el_in_blk1 = 3 ;
1954 num_fac_per_el1 = 17 ;
1955
1956 * num faces will be total face conn
1957 * num_faces_in_block will be number of faces (non repeated)
1958 * num_nodes_per_face will have total face connectivity
1959 */
1960 int num_faces2 = 0;
1961
1962 for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
1963 {
1964 EntityHandle polyh = *eit;
1965 int nfaces = 0;
1966 const EntityHandle* conn = NULL;
1967 ErrorCode rval = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
1968 num_faces2 += nfaces;
1969 }
1970
1971 int num_fac_per_el;
1972
1973 INS_ID( wname, "num_fac_per_el%d", element_block_index );
1974 if( nc_def_dim( ncFile, wname, (size_t)num_faces2, &num_fac_per_el ) != NC_NOERR )
1975 {
1976 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of faces per block " << block.id );
1977 }
1978
1979 /*
1980 int facconn1(num_fac_per_el1) ;
1981 facconn1:elem_type = "NFACED" ;
1982 int ebepecnt1(num_el_in_blk1) ;
1983 ebepecnt1:entity_type1 = "FACE" ;
1984 ebepecnt1:entity_type2 = "ELEM" ;
1985 */
1986
1987 // facconn
1988 INS_ID( wname, "facconn%d", element_block_index );
1989 int facconn;
1990 dims[0] = num_fac_per_el;
1991 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &facconn ) )
1992 {
1993 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create facconn array for block " << i + 1 );
1994 }
1995 std::string etype( "NFACED" );
1996 if( NC_NOERR != nc_put_att_text( ncFile, facconn, "elem_type", etype.length(), etype.c_str() ) )
1997 {
1998 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store elem type " << (int)block.element_type );
1999 }
2000
2001 // ebepecnt
2002 INS_ID( wname, "ebepecnt%d", element_block_index );
2003 int ebepecnt;
2004 dims[0] = num_el_in_blk;
2005 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
2006 {
2007 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
2008 }
2009 std::string etype1( "FACE" );
2010 if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
2011 {
2012 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
2013 }
2014 std::string etype2( "ELEM" );
2015 if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
2016 {
2017 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
2018 }
2019
2020 block.number_nodes_per_element = num_faces2; // connectivity for all polyhedra in block
2021 }
2022 }
2023
2024 /* Node set id array: */
2025
2026 int non_empty_nss = 0;
2027 // Need to go through nodesets to compute # nodesets, some might be empty
2028 std::vector< DirichletSetData >::iterator ns_it;
2029 for( ns_it = nodeset_data.begin(); ns_it != nodeset_data.end(); ++ns_it )
2030 {
2031 if( 0 != ( *ns_it ).number_nodes ) non_empty_nss++;
2032 }
2033
2034 int num_ns = -1;
2035 int ns_idstat = -1, ns_idarr = -1;
2036 if( non_empty_nss > 0 )
2037 {
2038 if( nc_def_dim( ncFile, "num_node_sets", (size_t)( non_empty_nss ), &num_ns ) != NC_NOERR )
2039 {
2040 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of node sets" );
2041 }
2042
2043 /* Node set id status array: */
2044
2045 if( NC_NOERR != nc_def_var( ncFile, "ns_status", NC_LONG, 1, &num_ns, &ns_idstat ) )
2046 {
2047 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets status array" );
2048 }
2049
2050 /* Node set id array: */
2051 if( NC_NOERR != nc_def_var( ncFile, "ns_prop1", NC_LONG, 1, &num_ns, &ns_idarr ) )
2052 {
2053 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets property array" );
2054 }
2055
2056 /* Store property name as attribute of property array variable */
2057 if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "name", strlen( "ID" ), "ID" ) )
2058 {
2059 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store node set property name ID" );
2060 }
2061
2062 // Now, define the arrays needed for each node set
2063
2064 int index = 0;
2065
2066 for( unsigned i = 0; i < nodeset_data.size(); i++ )
2067 {
2068 DirichletSetData node_set = nodeset_data[i];
2069
2070 if( 0 == node_set.number_nodes )
2071 {
2072 MB_SET_ERR_CONT( "WriteNCDF: empty nodeset " << node_set.id );
2073 continue;
2074 }
2075 index++;
2076
2077 int num_nod_ns = -1;
2078 INS_ID( wname, "num_nod_ns%d", index );
2079 if( nc_def_dim( ncFile, wname, (size_t)node_set.number_nodes, &num_nod_ns ) != NC_NOERR )
2080 {
2081 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for set " << node_set.id );
2082 }
2083
2084 /* Create variable array in which to store the node set node list */
2085 int node_ns = -1;
2086 INS_ID( wname, "node_ns%d", index );
2087 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_nod_ns, &node_ns ) )
2088 {
2089 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node set " << node_set.id << " node list" );
2090 }
2091
2092 // Create distribution factor array
2093 int fact_ns = -1;
2094 INS_ID( wname, "dist_fact_ns%d", index );
2095 if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 1, &num_nod_ns, &fact_ns ) )
2096 {
2097 MB_SET_ERR( MB_FAILURE,
2098 "WriteNCDF: failed to create node set " << node_set.id << " distribution factor list" );
2099 }
2100 }
2101 }
2102
2103 /* Side set id array: */
2104
2105 long non_empty_ss = 0;
2106 // Need to go through nodesets to compute # nodesets, some might be empty
2107 std::vector< NeumannSetData >::iterator ss_it;
2108 for( ss_it = sideset_data.begin(); ss_it != sideset_data.end(); ++ss_it )
2109 {
2110 if( 0 != ( *ss_it ).number_elements ) non_empty_ss++;
2111 }
2112
2113 if( non_empty_ss > 0 )
2114 {
2115 int num_ss = -1;
2116 if( nc_def_dim( ncFile, "num_side_sets", non_empty_ss, &num_ss ) != NC_NOERR )
2117 {
2118 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of side sets" );
2119 }
2120
2121 /* Side set id status array: */
2122 int ss_idstat = -1, ss_idarr = -1;
2123 if( NC_NOERR != nc_def_var( ncFile, "ss_status", NC_LONG, 1, &num_ss, &ss_idstat ) )
2124 {
2125 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set status" );
2126 }
2127
2128 /* Side set id array: */
2129 if( NC_NOERR != nc_def_var( ncFile, "ss_prop1", NC_LONG, 1, &num_ss, &ss_idarr ) )
2130 {
2131 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set property" );
2132 }
2133
2134 /* Store property name as attribute of property array variable */
2135 if( NC_NOERR != nc_put_att_text( ncFile, ss_idarr, "name", strlen( "ID" ), "ID" ) )
2136 {
2137 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store side set property name ID" );
2138 }
2139
2140 // Now, define the arrays needed for each side set
2141
2142 int index = 0;
2143 for( unsigned int i = 0; i < sideset_data.size(); i++ )
2144 {
2145 NeumannSetData side_set = sideset_data[i];
2146
2147 // Don't define an empty set
2148 if( 0 == side_set.number_elements ) continue;
2149
2150 index++;
2151
2152 int num_side_ss = -1;
2153 int elem_ss = -1, side_ss = -1;
2154 INS_ID( wname, "num_side_ss%d", index );
2155 if( nc_def_dim( ncFile, wname, (size_t)side_set.number_elements, &num_side_ss ) != NC_NOERR )
2156 {
2157 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of sides in side set " << side_set.id );
2158 }
2159
2160 INS_ID( wname, "elem_ss%d", index );
2161 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &elem_ss ) )
2162 {
2163 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element list for side set "
2164 << side_set.id ); /* Exit define mode and return */
2165 }
2166 INS_ID( wname, "side_ss%d", index );
2167 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &side_ss ) )
2168 {
2169 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create side list for side set "
2170 << side_set.id ); /* Exit define mode and return */
2171 }
2172
2173 // sideset distribution factors
2174 int num_df_ss = -1;
2175 INS_ID( wname, "num_df_ss%d", index );
2176 if( nc_def_dim( ncFile, wname, (size_t)side_set.ss_dist_factors.size(), &num_df_ss ) != NC_NOERR )
2177 {
2178 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dist factors in side set "
2179 << side_set.id ); /* Exit define mode and return */
2180 }
2181
2182 /* Create variable array in which to store the side set distribution factors */
2183
2184 int fact_ss = -1;
2185 INS_ID( wname, "dist_fact_ss%d", index );
2186 if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_df_ss, &fact_ss ) )
2187 {
2188 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create dist factors list for side set "
2189 << side_set.id ); /* Exit define mode and return */
2190 }
2191 }
2192 }
2193
2194 /* Node coordinate arrays: */
2195
2196 int coord, name_coord, dims[3];
2197 dims[0] = num_dim;
2198 dims[1] = num_nodes;
2199 if( NC_NOERR != nc_def_var( ncFile, "coord", NC_DOUBLE, 2, dims, &coord ) )
2200 {
2201 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define node coordinate array" );
2202 }
2203
2204 /* Coordinate names array */
2205
2206 dims[0] = num_dim;
2207 dims[1] = dim_str;
2208 if( NC_NOERR != nc_def_var( ncFile, "coor_names", NC_CHAR, 2, dims, &name_coord ) )
2209 {
2210 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define coordinate name array" );
2211 }
2212
2213 // Define genesis maps if required
2214
2215 if( write_maps )
2216 {
2217 // Element map
2218 int elem_map = -1, elem_map2 = -1, node_map = -1;
2219 if( NC_NOERR != nc_def_var( ncFile, "elem_map", NC_LONG, 1, &num_elem, &elem_map ) )
2220 {
2221 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element map array" ); /* Exit define mode and return */
2222 }
2223
2224 // Create the element number map
2225 if( NC_NOERR != nc_def_var( ncFile, "elem_num_map", NC_LONG, 1, &num_elem, &elem_map2 ) )
2226 {
2227 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element numbering map" ); /* Exit define mode
2228 and return */
2229 }
2230
2231 // Create node number map
2232 if( NC_NOERR != nc_def_var( ncFile, "node_num_map", NC_LONG, 1, &num_nodes, &node_map ) )
2233 {
2234 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node numbering map array" ); /* Exit define mode and
2235 return */
2236 }
2237 }
2238
2239 // Define qa records to be used
2240
2241 int num_qa_rec = mesh_info.qaRecords.size() / 4;
2242 int num_qa = -1;
2243
2244 if( nc_def_dim( ncFile, "num_qa_rec", (long)num_qa_rec, &num_qa ) != NC_NOERR )
2245 {
2246 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array size" );
2247 }
2248
2249 // Define qa array
2250 int qa_title;
2251 dims[0] = num_qa;
2252 dims[1] = dim_four;
2253 dims[2] = dim_str;
2254 if( NC_NOERR != nc_def_var( ncFile, "qa_records", NC_CHAR, 3, dims, &qa_title ) )
2255 {
2256 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array" );
2257 }
2258
2259 // Take it out of define mode
2260 if( NC_NOERR != nc_enddef( ncFile ) )
2261 {
2262 MB_SET_ERR( MB_FAILURE, "WriteNCDF: Trouble leaving define mode" );
2263 }
2264
2265 return MB_SUCCESS;
2266 }
2267
2268 ErrorCode WriteNCDF::open_file( const char* filename )
2269 {
2270 // Not a valid filename
2271 if( strlen( (const char*)filename ) == 0 )
2272 {
2273 MB_SET_ERR( MB_FAILURE, "Output Exodus filename not specified" );
2274 }
2275
2276 int fail = nc_create( filename, NC_CLOBBER, &ncFile );
2277
2278 // File couldn't be opened
2279 if( NC_NOERR != fail )
2280 {
2281 MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
2282 }
2283
2284 return MB_SUCCESS;
2285 }
2286
2287 ErrorCode WriteNCDF::get_sideset_elems( EntityHandle sideset,
2288 int current_sense,
2289 Range& forward_elems,
2290 Range& reverse_elems )
2291 {
2292 Range ss_elems, ss_meshsets;
2293
2294 // Get the sense tag; don't need to check return, might be an error if the tag
2295 // hasn't been created yet
2296 Tag sense_tag = 0;
2297 mdbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
2298
2299 // Get the entities in this set
2300 ErrorCode result = mdbImpl->get_entities_by_handle( sideset, ss_elems, true );
2301 if( MB_FAILURE == result ) return result;
2302
2303 // Now remove the meshsets into the ss_meshsets; first find the first meshset,
2304 Range::iterator range_iter = ss_elems.begin();
2305 while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
2306 ++range_iter;
2307
2308 // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
2309 if( range_iter != ss_elems.end() )
2310 {
2311 std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
2312 ss_elems.erase( range_iter, ss_elems.end() );
2313 }
2314
2315 // OK, for the elements, check the sense of this set and copy into the right range
2316 // (if the sense is 0, copy into both ranges)
2317
2318 // Need to step forward on list until we reach the right dimension
2319 Range::iterator dum_it = ss_elems.end();
2320 --dum_it;
2321 int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
2322 dum_it = ss_elems.begin();
2323 while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
2324 ++dum_it;
2325
2326 if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
2327 if( current_sense == -1 || current_sense == 0 )
2328 std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
2329
2330 // Now loop over the contained meshsets, getting the sense of those and calling this
2331 // function recursively
2332 for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
2333 {
2334 // First get the sense; if it's not there, by convention it's forward
2335 int this_sense;
2336 if( 0 == sense_tag || MB_FAILURE == mdbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
2337 this_sense = 1;
2338
2339 // Now get all the entities on this meshset, with the proper (possibly reversed) sense
2340 get_sideset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
2341 }
2342
2343 return result;
2344 }
2345
2346 } // namespace moab
1.9.1.