Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
WriteNCDF.cpp
Go to the documentation of this file.
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  MB_CHK_ERR( mdbImpl->get_connectivity( polg, conn, nnodes ) );
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  MB_CHK_ERR( mdbImpl->get_connectivity( block_data.elements, mesh_info.polyhedronFaces ) );
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 );
607  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
608 
609  if( 0x1 == node_marked )
610  {
611  nodeset_data.nodes.push_back( *iter );
612  if( 0 != has_dist_factors )
613  nodeset_data.node_dist_factors.push_back( dist_factor_vector[j] );
614  else
615  nodeset_data.node_dist_factors.push_back( 1.0 );
616  }
617  j++;
618  }
619 
620  nodeset_data.number_nodes = nodeset_data.nodes.size();
621  nodeset_info.push_back( nodeset_data );
622  }
623 
624  //------sidesets--------
625  vector_iter = sidesets.begin();
626  end_vector_iter = sidesets.end();
627 
628  for( ; vector_iter != end_vector_iter; ++vector_iter )
629  {
630  NeumannSetData sideset_data;
631 
632  // Get the sideset's id
633  if( mdbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
634 
635  sideset_data.id = id;
636  sideset_data.mesh_set_handle = *vector_iter;
637 
638  // Get the sides in two lists, one forward the other reverse; starts with forward sense
639  // by convention
640  Range forward_elems, reverse_elems;
641  if( get_sideset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
642 
643  ErrorCode result = get_valid_sides( forward_elems, mesh_info, 1, sideset_data );
644  MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
645  result = get_valid_sides( reverse_elems, mesh_info, -1, sideset_data );
646  MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
647 
648  sideset_data.number_elements = sideset_data.elements.size();
649  sideset_info.push_back( sideset_data );
650  }
651 
652  return MB_SUCCESS;
653 }
654 
655 ErrorCode WriteNCDF::get_valid_sides( Range& elems,
656  ExodusMeshInfo& /*mesh_info*/,
657  const int sense,
658  NeumannSetData& sideset_data )
659 {
660  // This is where we see if underlying element of side set element is included in output
661 
662  // Get the sideset-based info for distribution factors
663  const double* dist_factor_vector = 0;
664  int dist_factor_size = 0;
665 
666  // Initialize dist_fac_iter to get rid of compiler warning
667  const double* dist_fac_iter = 0;
668  const void* ptr = 0;
669  bool has_dist_factors = false;
670  if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( sideset_data.mesh_set_handle ), 1, &ptr, &dist_factor_size ) ==
671  MB_SUCCESS &&
672  dist_factor_size )
673  {
674  has_dist_factors = true;
675  dist_factor_vector = reinterpret_cast< const double* >( ptr );
676  dist_fac_iter = dist_factor_vector;
677  dist_factor_size /= sizeof( double );
678  }
679 
680  unsigned char element_marked = 0;
681  ErrorCode result;
682  for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
683  {
684  // Should insert here if "side" is a quad/tri on a quad/tri mesh
685  result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );
686  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
687 
688  if( 0x1 == element_marked )
689  {
690  sideset_data.elements.push_back( *iter );
691 
692  // TJT TODO: the sense should really be # edges + 1or2
693  sideset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
694  }
695  else
696  { // Then "side" is probably a quad/tri on a hex/tet mesh
697  std::vector< EntityHandle > parents;
698  int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
699 
700  // Get the adjacent parent element of "side"
701  if( mdbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
702  {
703 #if 0
704  // This is not treated as an error, print warning messages for
705  // debugging only
706  fprintf(stderr, "[Warning]: Couldn't get adjacencies for sideset.\n");
707 #endif
708  }
709 
710  if( !parents.empty() )
711  {
712  // Make sure the adjacent parent element will be output
713  for( unsigned int k = 0; k < parents.size(); k++ )
714  {
715  result = mdbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );
716  MB_CHK_SET_ERR( result, "Couldn't get mark data" );
717 
718  int side_no, this_sense, this_offset;
719  if( 0x1 == element_marked &&
720  mdbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
721  this_sense == sense )
722  {
723  sideset_data.elements.push_back( parents[k] );
724  sideset_data.side_numbers.push_back( side_no + 1 );
725  break;
726  }
727  }
728  }
729  else
730  {
731 #if 0
732  // This is not treated as an error, print warning messages for
733  // debugging only
734  fprintf(stderr, "[Warning]: No parent element exists for element in sideset %i\n", sideset_data.id);
735 #endif
736  }
737  }
738 
739  if( sideset_data.elements.size() != 0 )
740  {
741  // Distribution factors
742  int num_nodes = CN::VerticesPerEntity( TYPE_FROM_HANDLE( *iter ) );
743  // put some dummy dist factors for polygons; why do we need them?
744  if( TYPE_FROM_HANDLE( *iter ) == MBPOLYGON ) num_nodes = 1; // dummy
745  if( has_dist_factors )
746  {
747  std::copy( dist_fac_iter, dist_fac_iter + num_nodes,
748  std::back_inserter( sideset_data.ss_dist_factors ) );
749  dist_fac_iter += num_nodes;
750  }
751  else
752  {
753  for( int j = 0; j < num_nodes; j++ )
754  sideset_data.ss_dist_factors.push_back( 1.0 );
755  }
756  }
757  }
758 
759  return MB_SUCCESS;
760 }
761 
762 ErrorCode WriteNCDF::write_qa_records( std::vector< std::string >& qa_record_list )
763 {
764  int i = 0;
765 
766  for( std::vector< std::string >::iterator string_it = qa_record_list.begin(); string_it != qa_record_list.end(); )
767  {
768  for( int j = 0; j < 4; j++ )
769  write_qa_string( ( *string_it++ ).c_str(), i, j );
770  i++;
771  }
772 
773  return MB_SUCCESS;
774 }
775 
776 ErrorCode WriteNCDF::write_qa_string( const char* string, int record_number, int record_position )
777 {
778  // Get the variable id in the exodus file
779 
780  std::vector< int > dims;
781  int temp_var = -1;
782  GET_VAR( "qa_records", temp_var, dims );
783  if( -1 == temp_var )
784  {
785  MB_SET_ERR( MB_FAILURE, "WriteNCDF:: Problem getting qa record variable" );
786  }
787  size_t count[3], start[3];
788 
789  // Write out the record
790  start[0] = record_number;
791  start[1] = record_position;
792  start[2] = 0;
793 
794  count[0] = 1;
795  count[1] = 1;
796  count[2] = (long)strlen( string ) + 1;
797  int fail = nc_put_vara_text( ncFile, temp_var, start, count, string );
798  if( NC_NOERR != fail )
799  {
800  MB_SET_ERR( MB_FAILURE, "Failed to position qa string variable" );
801  }
802 
803  return MB_SUCCESS;
804 }
805 
806 ErrorCode WriteNCDF::write_nodes( int num_nodes, Range& nodes, int dimension )
807 {
808  // Write coordinates names
809  int nc_var = -1;
810  std::vector< int > dims;
811  GET_VAR( "coor_names", nc_var, dims );
812  if( -1 == nc_var )
813  {
814  MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate name variable" );
815  }
816 
817  size_t start[2] = { 0, 0 }, count[2] = { 1, ExoIIInterface::MAX_STR_LENGTH };
818  char dum_str[ExoIIInterface::MAX_STR_LENGTH];
819  strcpy( dum_str, "x" );
820  int fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
821  if( NC_NOERR != fail )
822  {
823  MB_SET_ERR( MB_FAILURE, "Trouble adding x coordinate name; netcdf message: " << nc_strerror( fail ) );
824  }
825 
826  start[0] = 1;
827  strcpy( dum_str, "y" );
828  fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
829  if( NC_NOERR != fail )
830  {
831  MB_SET_ERR( MB_FAILURE, "Trouble adding y coordinate name; netcdf message: " << nc_strerror( fail ) );
832  }
833 
834  start[0] = 2;
835  strcpy( dum_str, "z" );
836  fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
837  if( NC_NOERR != fail )
838  {
839  MB_SET_ERR( MB_FAILURE, "Trouble adding z coordinate name; netcdf message: " << nc_strerror( fail ) );
840  }
841 
842  // See if should transform coordinates
843  ErrorCode result;
844  Tag trans_tag;
845  result = mdbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
846  bool transform_needed = true;
847  if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
848 
849  int num_coords_to_fill = transform_needed ? 3 : dimension;
850 
851  std::vector< double* > coord_arrays( 3 );
852  coord_arrays[0] = new double[num_nodes];
853  coord_arrays[1] = new double[num_nodes];
854  coord_arrays[2] = NULL;
855 
856  if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
857 
858  result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 1, coord_arrays );
859  if( result != MB_SUCCESS )
860  {
861  delete[] coord_arrays[0];
862  delete[] coord_arrays[1];
863  if( coord_arrays[2] ) delete[] coord_arrays[2];
864  return result;
865  }
866 
867  if( transform_needed )
868  {
869  double trans_matrix[16];
870  const EntityHandle mesh = 0;
871  result = mdbImpl->tag_get_data( trans_tag, &mesh, 0, trans_matrix );
872  MB_CHK_SET_ERR( result, "Couldn't get transform data" );
873 
874  for( int i = 0; i < num_nodes; i++ )
875  {
876  double vec1[3];
877  double vec2[3];
878 
879  vec2[0] = coord_arrays[0][i];
880  vec2[1] = coord_arrays[1][i];
881  vec2[2] = coord_arrays[2][i];
882 
883  for( int row = 0; row < 3; row++ )
884  {
885  vec1[row] = 0.0;
886  for( int col = 0; col < 3; col++ )
887  vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
888  }
889 
890  coord_arrays[0][i] = vec1[0];
891  coord_arrays[1][i] = vec1[1];
892  coord_arrays[2][i] = vec1[2];
893  }
894  }
895 
896  // Write the nodes
897  nc_var = -1;
898  GET_VAR( "coord", nc_var, dims );
899  if( -1 == nc_var )
900  {
901  MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate variable" );
902  }
903  start[0] = 0;
904  count[1] = num_nodes;
905  fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[0][0] ) );
906  if( NC_NOERR != fail )
907  {
908  MB_SET_ERR( MB_FAILURE, "Trouble writing x coordinate" );
909  }
910 
911  start[0] = 1;
912  fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[1][0] ) );
913  if( NC_NOERR != fail )
914  {
915  MB_SET_ERR( MB_FAILURE, "Trouble writing y coordinate" );
916  }
917 
918  start[0] = 2;
919  fail = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[2][0] ) );
920  if( NC_NOERR != fail )
921  {
922  MB_SET_ERR( MB_FAILURE, "Trouble writing z coordinate" );
923  }
924 
925  delete[] coord_arrays[0];
926  delete[] coord_arrays[1];
927  if( coord_arrays[2] ) delete[] coord_arrays[2];
928 
929  return MB_SUCCESS;
930 }
931 
932 ErrorCode WriteNCDF::write_poly_faces( ExodusMeshInfo& mesh_info )
933 {
934  // write all polygons that are not in another element block;
935  // usually they are nowhere else, but be sure, write in this block only ones that are not in the
936  // other blocks
937  Range pfaces = mesh_info.polyhedronFaces;
938 
939  /*
940  * int fbconn1(num_nod_per_fa1) ;
941  fbconn1:elem_type = "nsided" ;
942  int fbepecnt1(num_fa_in_blk1) ;
943  fbepecnt1:entity_type1 = "NODE" ;
944  fbepecnt1:entity_type2 = "FACE" ;
945  */
946  if( pfaces.empty() ) return MB_SUCCESS;
947  char wname[CHAR_STR_LEN];
948  int nc_var = -1;
949  std::vector< int > dims;
950 
951  // write one for each element block, to make paraview and visit happy
952  int num_faces_in_block = (int)pfaces.size();
953  for( unsigned int bl = 0; bl < mesh_info.num_polyhedra_blocks; bl++ )
954  {
955  INS_ID( wname, "fbconn%u", bl + 1 ); // it is the first block
956  GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
957 
958  INS_ID( wname, "num_nod_per_fa%u", bl + 1 );
959  int ncdim, num_nod_per_face;
960  GET_DIM( ncdim, wname, num_nod_per_face );
961  int* connectivity = new int[num_nod_per_face];
962  int ixcon = 0, j = 0;
963  std::vector< int > fbepe( num_faces_in_block ); // fbepecnt1
964  for( Range::iterator eit = pfaces.begin(); eit != pfaces.end(); eit++ )
965  {
966  EntityHandle polyg = *eit;
967  int nnodes = 0;
968  const EntityHandle* conn = NULL;
969  MB_CHK_ERR( mdbImpl->get_connectivity( polyg, conn, nnodes ) );
970  for( int k = 0; k < nnodes; k++ )
971  connectivity[ixcon++] = conn[k];
972  fbepe[j++] = nnodes;
973  }
974  size_t start[1] = { 0 }, count[1] = { 0 };
975  count[0] = ixcon;
976  int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
977  if( NC_NOERR != fail )
978  {
979  delete[] connectivity;
980  MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
981  }
982 
983  INS_ID( wname, "fbepecnt%u", bl + 1 );
984  GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
985  count[0] = num_faces_in_block;
986 
987  fail = nc_put_vara_int( ncFile, nc_var, start, count, &fbepe[0] );
988  if( NC_NOERR != fail )
989  {
990  delete[] connectivity;
991  MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
992  }
993 
994  int id = bl + 1;
995  if( write_exodus_integer_variable( "fa_prop1", &id, bl, 1 ) != MB_SUCCESS )
996  {
997  MB_SET_ERR_CONT( "Problem writing element block id " << id );
998  }
999 
1000  int status = 1;
1001  if( write_exodus_integer_variable( "fa_status", &status, bl, 1 ) != MB_SUCCESS )
1002  {
1003  MB_SET_ERR( MB_FAILURE, "Problem writing face block status" );
1004  }
1005 
1006  delete[] connectivity;
1007  if( 0 == repeat_face_blocks ) break; // do not repeat face blocks
1008  }
1009 
1010  return MB_SUCCESS;
1011 }
1012 ErrorCode WriteNCDF::write_header( ExodusMeshInfo& mesh_info,
1013  std::vector< MaterialSetData >& block_info,
1014  std::vector< NeumannSetData >& sideset_info,
1015  std::vector< DirichletSetData >& nodeset_info,
1016  const char* filename )
1017 {
1018  // Get the date and time
1019  char time[TIME_STR_LEN];
1020  char date[TIME_STR_LEN];
1021  time_and_date( time, date );
1022 
1023  std::string title_string = "MOAB";
1024  title_string.append( "(" );
1025  title_string.append( filename );
1026  title_string.append( "): " );
1027  title_string.append( date );
1028  title_string.append( ": " );
1029  title_string.append( "time " );
1030 
1031  if( title_string.length() > ExoIIInterface::MAX_LINE_LENGTH )
1032  title_string.resize( ExoIIInterface::MAX_LINE_LENGTH );
1033 
1034  // Initialize the exodus file
1035 
1036  int result = initialize_exodus_file( mesh_info, block_info, sideset_info, nodeset_info, title_string.c_str() );
1037 
1038  if( result == MB_FAILURE ) return MB_FAILURE;
1039 
1040  return MB_SUCCESS;
1041 }
1042 
1043 ErrorCode WriteNCDF::write_elementblocks( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data )
1044 {
1045  unsigned int i;
1046  int block_index = 0; // Index into block list, may != 1 if there are inactive blocks
1047  int exodus_id = 1;
1048 
1049  for( i = 0; i < block_data.size(); i++ )
1050  {
1051  MaterialSetData& block = block_data[i];
1052 
1053  unsigned int num_nodes_per_elem = block.number_nodes_per_element;
1054 
1055  // Write out the id for the block
1056 
1057  int id = block.id;
1058  int num_values = 1;
1059 
1060  if( write_exodus_integer_variable( "eb_prop1", &id, block_index, num_values ) != MB_SUCCESS )
1061  {
1062  MB_SET_ERR_CONT( "Problem writing element block id " << id );
1063  }
1064 
1065  // Write out the block status
1066 
1067  int status = 1;
1068  if( 0 == block.number_elements )
1069  {
1070  MB_SET_ERR( MB_FAILURE, "No elements in block " << id );
1071  }
1072 
1073  if( write_exodus_integer_variable( "eb_status", &status, block_index, num_values ) != MB_SUCCESS )
1074  {
1075  MB_SET_ERR( MB_FAILURE, "Problem writing element block status" );
1076  }
1077 
1078  //
1079  // Map the connectivity to the new nodes
1080  const unsigned int num_elem = block.number_elements;
1081  unsigned int num_nodes = num_nodes_per_elem * num_elem;
1082  if( EXOII_POLYGON == block.element_type || EXOII_POLYHEDRON == block.element_type )
1083  {
1084  num_nodes = num_nodes_per_elem;
1085  }
1086  int* connectivity = new int[num_nodes];
1087 
1088  ErrorCode result = MB_SUCCESS;
1089  if( block.element_type != EXOII_POLYHEDRON )
1090  mWriteIface->get_element_connect( num_elem, num_nodes_per_elem, mGlobalIdTag, block.elements, mGlobalIdTag,
1091  exodus_id, connectivity );
1092  if( result != MB_SUCCESS )
1093  {
1094  delete[] connectivity;
1095  MB_SET_ERR( result, "Couldn't get element array to write from" );
1096  }
1097 
1098  // If necessary, convert from EXODUS to CN node order
1099  const EntityType elem_type = ExoIIUtil::ExoIIElementMBEntity[block.element_type];
1100  assert( block.elements.all_of_type( elem_type ) );
1101  const int* reorder = 0;
1102  if( block.element_type != EXOII_POLYHEDRON && block.element_type != EXOII_POLYGON )
1103  reorder = exodus_elem_order_map[elem_type][block.number_nodes_per_element];
1104  if( reorder )
1105  WriteUtilIface::reorder( reorder, connectivity, block.number_elements, block.number_nodes_per_element );
1106 
1107  char wname[CHAR_STR_LEN];
1108  int nc_var = -1;
1109  std::vector< int > dims;
1110  if( block.element_type != EXOII_POLYHEDRON )
1111  {
1112  exodus_id += num_elem;
1113  INS_ID( wname, "connect%u", i + 1 );
1114 
1115  GET_VAR( wname, nc_var, dims );
1116  if( -1 == nc_var )
1117  {
1118  delete[] connectivity;
1119  MB_SET_ERR( MB_FAILURE, "Couldn't get connectivity variable" );
1120  }
1121  }
1122 
1123  if( EXOII_POLYGON == block.element_type )
1124  {
1125  size_t start[1] = { 0 }, count[1] = { num_nodes_per_elem };
1126  int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1127  if( NC_NOERR != fail )
1128  {
1129  delete[] connectivity;
1130  MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
1131  }
1132  // now put also number ebepecnt1
1133  INS_ID( wname, "ebepecnt%u", i + 1 );
1134  GET_VAR( wname, nc_var, dims );
1135  count[0] = block.number_elements;
1136  start[0] = 0;
1137  // reuse connectivity array, to not allocate another one
1138  int j = 0;
1139  for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); j++, eit++ )
1140  {
1141  EntityHandle polg = *eit;
1142  int nnodes = 0;
1143  const EntityHandle* conn = NULL;
1144  MB_CHK_ERR( mdbImpl->get_connectivity( polg, conn, nnodes ) );
1145  connectivity[j] = nnodes;
1146  }
1147  fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1148  if( NC_NOERR != fail )
1149  {
1150  delete[] connectivity;
1151  MB_SET_ERR( MB_FAILURE, "Couldn't write ebepecnt variable" );
1152  }
1153  }
1154  else if( block.element_type != EXOII_POLYHEDRON )
1155  {
1156  size_t start[2] = { 0, 0 }, count[2] = { num_elem, num_nodes_per_elem };
1157  int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1158  if( NC_NOERR != fail )
1159  {
1160  delete[] connectivity;
1161  MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
1162  }
1163  }
1164  else // if (block.element_type == EXOII_POLYHEDRON)
1165  {
1166  /* write a lot of stuff // faconn
1167  num_fa_in_blk1 = 15 ;
1168  num_nod_per_fa1 = 58 ;
1169  num_el_in_blk1 = 3 ;
1170  num_fac_per_el1 = 17 ;
1171  int fbconn1(num_nod_per_fa1) ;
1172  fbconn1:elem_type = "nsided" ;
1173  int fbepecnt1(num_fa_in_blk1) ;
1174  fbepecnt1:entity_type1 = "NODE" ;
1175  fbepecnt1:entity_type2 = "FACE" ;
1176  int facconn1(num_fac_per_el1) ;
1177  facconn1:elem_type = "NFACED" ;
1178  int ebepecnt1(num_el_in_blk1) ;
1179  ebepecnt1:entity_type1 = "FACE" ;
1180  ebepecnt1:entity_type2 = "ELEM" ;
1181 
1182  */
1183  Range& block_faces = mesh_info.polyhedronFaces;
1184  // ErrorCode rval = mdbImpl->get_connectivity(block.elements, block_faces);
1185  // MB_CHK_ERR(rval);
1186 
1187  // reuse now connectivity for facconn1
1188  INS_ID( wname, "facconn%u", i + 1 );
1189  GET_VAR( wname, nc_var, dims ); // fbconn# variable, 1 dimensional
1190 
1191  std::vector< int > ebepe( block.elements.size() ); // ebepecnt1
1192  int ixcon = 0, j = 0;
1193  size_t start[1] = { 0 }, count[1] = { 0 };
1194 
1195  for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
1196  {
1197  EntityHandle polyh = *eit;
1198  int nfaces = 0;
1199  const EntityHandle* conn = NULL;
1200  MB_CHK_ERR( mdbImpl->get_connectivity( polyh, conn, nfaces ) );
1201  for( int k = 0; k < nfaces; k++ )
1202  {
1203  int index = block_faces.index( conn[k] );
1204  if( index == -1 ) MB_SET_ERR( MB_FAILURE, "Couldn't find face in polyhedron" );
1205  connectivity[ixcon++] = index + 1;
1206  }
1207  ebepe[j++] = nfaces;
1208  // num_faces+=nfaces;
1209  }
1210  count[0] = ixcon; // facconn1
1211  int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
1212  if( NC_NOERR != fail )
1213  {
1214  delete[] connectivity;
1215  MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
1216  }
1217 
1218  INS_ID( wname, "ebepecnt%u", i + 1 );
1219  GET_VAR( wname, nc_var, dims ); // ebepecnt# variable, 1 dimensional
1220  count[0] = block.elements.size();
1221 
1222  fail = nc_put_vara_int( ncFile, nc_var, start, count, &ebepe[0] );
1223  if( NC_NOERR != fail )
1224  {
1225  delete[] connectivity;
1226  MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
1227  }
1228  }
1229  block_index++;
1230  delete[] connectivity;
1231  }
1232 
1233  return MB_SUCCESS;
1234 }
1235 
1236 ErrorCode WriteNCDF::write_global_node_order_map( int num_nodes, Range& nodes )
1237 {
1238  // Note: this routine bypasses the standard exodusII interface for efficiency!
1239 
1240  // Node order map
1241  int* map = new int[num_nodes];
1242 
1243  // For now, output a dummy map!
1244 
1245  Range::iterator range_iter, end_iter;
1246  range_iter = nodes.begin();
1247  end_iter = nodes.end();
1248 
1249  int i = 0;
1250 
1251  for( ; range_iter != end_iter; ++range_iter )
1252  {
1253  // TODO -- do we really want to cast this to an int?
1254  map[i++] = (int)ID_FROM_HANDLE( *range_iter );
1255  }
1256 
1257  // Output array and cleanup
1258 
1259  int error = write_exodus_integer_variable( "node_num_map", map, 0, num_nodes );
1260 
1261  if( map ) delete[] map;
1262 
1263  if( error < 0 )
1264  {
1265  MB_SET_ERR( MB_FAILURE, "Failed writing global node order map" );
1266  }
1267 
1268  return MB_SUCCESS;
1269 }
1270 
1271 ErrorCode WriteNCDF::write_global_element_order_map( int num_elements )
1272 {
1273  // Allocate map array
1274  int* map = new int[num_elements];
1275 
1276  // Many Sandia codes assume this map is unique, and CUBIT does not currently
1277  // have unique ids for all elements. Therefore, to make sure nothing crashes,
1278  // insert a dummy map...
1279 
1280  for( int i = 0; i < num_elements; i++ )
1281  map[i] = i + 1;
1282 
1283  // Output array and cleanup
1284 
1285  int error = write_exodus_integer_variable( "elem_num_map", map, 0, num_elements );
1286 
1287  if( map ) delete[] map;
1288 
1289  if( error < 0 )
1290  {
1291  MB_SET_ERR( MB_FAILURE, "Failed writing global element order map" );
1292  }
1293 
1294  return MB_SUCCESS;
1295 }
1296 
1297 ErrorCode WriteNCDF::write_element_order_map( int num_elements )
1298 {
1299  // Note: this routine bypasses the standard exodusII interface for efficiency!
1300 
1301  // Element order map
1302  int* map = new int[num_elements];
1303 
1304  // For now, output a dummy map!
1305 
1306  for( int i = 0; i < num_elements; i++ )
1307  {
1308  map[i] = i + 1;
1309  }
1310 
1311  // Output array and cleanup
1312 
1313  int error = write_exodus_integer_variable( "elem_map", map, 0, num_elements );
1314 
1315  if( map ) delete[] map;
1316 
1317  if( error < 0 )
1318  {
1319  MB_SET_ERR( MB_FAILURE, "Failed writing element map" );
1320  }
1321 
1322  return MB_SUCCESS;
1323 }
1324 
1325 ErrorCode WriteNCDF::write_exodus_integer_variable( const char* variable_name,
1326  int* variable_array,
1327  int start_position,
1328  int number_values )
1329 {
1330  // Note: this routine bypasses the standard exodusII interface for efficiency!
1331 
1332  // Write directly to netcdf interface for efficiency
1333 
1334  // Get the variable id of the element map
1335  int nc_var = -1;
1336  std::vector< int > dims;
1337  GET_VAR( variable_name, nc_var, dims );
1338  if( -1 == nc_var )
1339  {
1340  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate variable " << variable_name << " in file" );
1341  }
1342  // This contortion is necessary because netCDF is expecting nclongs;
1343  // fortunately it's necessary only when ints and nclongs aren't the same size
1344 
1345  size_t start[1], count[1];
1346  start[0] = start_position;
1347  count[0] = number_values;
1348 
1349  int fail = NC_NOERR;
1350  if( sizeof( int ) == sizeof( long ) )
1351  {
1352  fail = nc_put_vara_int( ncFile, nc_var, start, count, variable_array );
1353  }
1354  else
1355  {
1356  long* lptr = new long[number_values];
1357  for( int jj = 0; jj < number_values; jj++ )
1358  lptr[jj] = variable_array[jj];
1359  fail = nc_put_vara_long( ncFile, nc_var, start, count, lptr );
1360  delete[] lptr;
1361  }
1362 
1363  if( NC_NOERR != fail )
1364  {
1365  MB_SET_ERR( MB_FAILURE, "Failed to store variable " << variable_name );
1366  }
1367 
1368  return MB_SUCCESS;
1369 }
1370 
1371 ErrorCode WriteNCDF::write_BCs( std::vector< NeumannSetData >& sidesets, std::vector< DirichletSetData >& nodesets )
1372 {
1373  unsigned int i, j;
1374  int id;
1375  int ns_index = -1;
1376 
1377  for( std::vector< DirichletSetData >::iterator ns_it = nodesets.begin(); ns_it != nodesets.end(); ++ns_it )
1378  {
1379  // Get number of nodes in set
1380  int number_nodes = ( *ns_it ).number_nodes;
1381  if( 0 == number_nodes ) continue;
1382 
1383  // If we're here, we have a non-empty nodeset; increment the index
1384  ns_index++;
1385 
1386  // Get the node set id
1387  id = ( *ns_it ).id;
1388 
1389  // Build new array to old exodus ids
1390  int* exodus_id_array = new int[number_nodes];
1391  double* dist_factor_array = new double[number_nodes];
1392 
1393  std::vector< EntityHandle >::iterator begin_iter, end_iter;
1394  std::vector< double >::iterator other_iter;
1395  begin_iter = ( *ns_it ).nodes.begin();
1396  end_iter = ( *ns_it ).nodes.end();
1397  other_iter = ( *ns_it ).node_dist_factors.begin();
1398 
1399  j = 0;
1400  int exodus_id;
1401  ErrorCode result;
1402  // Fill up node array and dist. factor array at the same time
1403  for( ; begin_iter != end_iter; ++begin_iter )
1404  {
1405  result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );
1406  MB_CHK_SET_ERR( result, "Problem getting id tag data" );
1407 
1408  exodus_id_array[j] = exodus_id;
1409  dist_factor_array[j] = *( other_iter );
1410  ++other_iter;
1411  j++;
1412  }
1413 
1414  // Write out the id for the nodeset
1415 
1416  int num_values = 1;
1417 
1418  result = write_exodus_integer_variable( "ns_prop1", &id, ns_index, num_values );
1419  MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
1420 
1421  // Write out the nodeset status
1422 
1423  int status = 1;
1424  if( !number_nodes ) status = 0;
1425 
1426  result = write_exodus_integer_variable( "ns_status", &status, ns_index, num_values );
1427  MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set status", MB_FAILURE );
1428 
1429  // Write it out
1430  char wname[CHAR_STR_LEN];
1431  int nc_var = -1;
1432  std::vector< int > dims;
1433  INS_ID( wname, "node_ns%d", ns_index + 1 );
1434  GET_VAR( wname, nc_var, dims );
1435  if( -1 == nc_var )
1436  {
1437  MB_SET_ERR( MB_FAILURE, "Failed to get node_ns variable" );
1438  }
1439 
1440  size_t start = 0, count = number_nodes;
1441  int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, exodus_id_array );
1442  if( NC_NOERR != fail )
1443  {
1444  MB_SET_ERR( MB_FAILURE, "Failed writing exodus id array" );
1445  }
1446 
1447  // Write out nodeset distribution factors
1448  INS_ID( wname, "dist_fact_ns%d", ns_index + 1 );
1449  nc_var = -1;
1450  GET_VAR( wname, nc_var, dims );
1451  if( -1 == nc_var )
1452  {
1453  MB_SET_ERR( MB_FAILURE, "Failed to get dist_fact variable" );
1454  }
1455  fail = nc_put_vara_double( ncFile, nc_var, &start, &count, dist_factor_array );
1456  if( NC_NOERR != fail )
1457  {
1458  MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" );
1459  }
1460 
1461  delete[] dist_factor_array;
1462  delete[] exodus_id_array;
1463  }
1464 
1465  // Now do sidesets
1466  int ss_index = 0; // Index of sideset - not the same as 'i' because
1467  // only writing non-empty side sets
1468  for( i = 0; i < sidesets.size(); i++ )
1469  {
1470  NeumannSetData sideset_data = sidesets[i];
1471 
1472  // Get the side set id
1473  int side_set_id = sideset_data.id;
1474 
1475  // Get number of elements in set
1476  int number_elements = sideset_data.number_elements;
1477  if( 0 == number_elements ) continue;
1478 
1479  // Build new array to old exodus ids
1480  int* output_element_ids = new int[number_elements];
1481  int* output_element_side_numbers = new int[number_elements];
1482 
1483  std::vector< EntityHandle >::iterator begin_iter, end_iter;
1484  begin_iter = sideset_data.elements.begin();
1485  end_iter = sideset_data.elements.end();
1486  std::vector< int >::iterator side_iter = sideset_data.side_numbers.begin();
1487 
1488  // Get the tag handle
1489  j = 0;
1490  int exodus_id;
1491 
1492  // For each "side"
1493  for( ; begin_iter != end_iter; ++begin_iter, ++side_iter )
1494  {
1495  ErrorCode result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );
1496  MB_CHK_SET_ERR( result, "Problem getting exodus id for sideset element "
1497  << (long unsigned int)ID_FROM_HANDLE( *begin_iter ) );
1498 
1499  output_element_ids[j] = exodus_id;
1500  output_element_side_numbers[j++] = *side_iter;
1501  }
1502 
1503  if( 0 != number_elements )
1504  {
1505  // Write out the id for the nodeset
1506 
1507  int num_values = 1;
1508 
1509  // ss_prop1[ss_index] = side_set_id
1510  ErrorCode result = write_exodus_integer_variable( "ss_prop1", &side_set_id, ss_index, num_values );
1511  MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
1512 
1513  // FIXME : Something seems wrong here. The we are within a block
1514  // started with if (0 != number_elements), so this condition is always
1515  // false. This code seems to indicate that we want to write all
1516  // sidesets, not just empty ones. But the code both here and in
1517  // initialize_exodus_file() skip empty side sets.
1518  // - j.k. 2007-03-09
1519  int status = 1;
1520  if( 0 == number_elements ) status = 0;
1521 
1522  // ss_status[ss_index] = status
1523  result = write_exodus_integer_variable( "ss_status", &status, ss_index, num_values );
1524  MB_CHK_SET_ERR_RET_VAL( result, "Problem writing side set status", MB_FAILURE );
1525 
1526  // Increment ss_index now because we want a) we need to
1527  // increment it somewhere within the if (0 != number_elements)
1528  // block and b) the above calls need a zero-based index
1529  // while the following use a one-based index.
1530  ++ss_index;
1531 
1532  char wname[CHAR_STR_LEN];
1533  int nc_var;
1534  std::vector< int > dims;
1535  INS_ID( wname, "elem_ss%d", ss_index );
1536  GET_VAR( wname, nc_var, dims );
1537  if( -1 == nc_var )
1538  {
1539  MB_SET_ERR( MB_FAILURE, "Failed to get elem_ss variable" );
1540  }
1541  size_t start = 0, count = number_elements;
1542  int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_ids );
1543  if( NC_NOERR != fail )
1544  {
1545  MB_SET_ERR( MB_FAILURE, "Failed writing sideset element array" );
1546  }
1547 
1548  INS_ID( wname, "side_ss%d", ss_index );
1549  nc_var = -1;
1550  GET_VAR( wname, nc_var, dims );
1551  if( -1 == nc_var )
1552  {
1553  MB_SET_ERR( MB_FAILURE, "Failed to get side_ss variable" );
1554  }
1555  fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_side_numbers );
1556  if( NC_NOERR != fail )
1557  {
1558  MB_SET_ERR( MB_FAILURE, "Failed writing sideset side array" );
1559  }
1560 
1561  INS_ID( wname, "dist_fact_ss%d", ss_index );
1562  nc_var = -1;
1563  GET_VAR( wname, nc_var, dims );
1564  if( -1 == nc_var )
1565  {
1566  MB_SET_ERR( MB_FAILURE, "Failed to get sideset dist factors variable" );
1567  }
1568  count = sideset_data.ss_dist_factors.size();
1569  fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &( sideset_data.ss_dist_factors[0] ) );
1570  if( NC_NOERR != fail )
1571  {
1572  MB_SET_ERR( MB_FAILURE, "Failed writing sideset dist factors array" );
1573  }
1574  }
1575 
1576  delete[] output_element_ids;
1577  delete[] output_element_side_numbers;
1578  }
1579 
1580  return MB_SUCCESS;
1581 }
1582 
1583 ErrorCode WriteNCDF::initialize_exodus_file( ExodusMeshInfo& mesh_info,
1584  std::vector< MaterialSetData >& block_data,
1585  std::vector< NeumannSetData >& sideset_data,
1586  std::vector< DirichletSetData >& nodeset_data,
1587  const char* title_string,
1588  bool write_maps,
1589  bool /* write_sideset_distribution_factors */ )
1590 {
1591  // This routine takes the place of the exodusii routine ex_put_init,
1592  // and additionally pre-defines variables such as qa, element blocks,
1593  // sidesets and nodesets in a single pass.
1594  //
1595  // This is necessary because of the way exodusII works. Each time the
1596  // netcdf routine endef is called to take the file out of define mode,
1597  // the entire file is copied before the new information is added.
1598  //
1599  // With very large files, this is simply not workable. This routine takes
1600  // the definition portions of all applicable exodus routines and puts them
1601  // in a single definition, preventing repeated copying of the file.
1602  //
1603  // Most of the code is copied directly from the applicable exodus routine,
1604  // and thus this routine may not seem very consistent in usage or variable
1605  // naming, etc.
1606 
1607  // Perform the initializations
1608 
1609  int element_block_index;
1610 
1611  // Inquire on defined string dimension and general dimension for qa
1612 
1613  int dim_str, dim_four, dim_line, dim_time;
1614  if( nc_def_dim( ncFile, "len_string", ExoIIInterface::MAX_STR_LENGTH, &dim_str ) != NC_NOERR )
1615  {
1616  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get string length in file" );
1617  }
1618 
1619  if( nc_def_dim( ncFile, "len_line", ExoIIInterface::MAX_STR_LENGTH, &dim_line ) != NC_NOERR )
1620  {
1621  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get line length in file" );
1622  }
1623 
1624  if( nc_def_dim( ncFile, "four", 4, &dim_four ) != NC_NOERR )
1625  {
1626  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate four in file" );
1627  }
1628 
1629  if( nc_def_dim( ncFile, "time_step", 1, &dim_time ) != NC_NOERR )
1630  {
1631  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate time step in file" );
1632  }
1633  // some whole_time dummy :(
1634  int dtime;
1635  if( NC_NOERR != nc_def_var( ncFile, "time_whole", NC_DOUBLE, 1, &dim_time, &dtime ) )
1636  {
1637  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define time whole array" );
1638  }
1639 
1640  /* Put file into define mode */
1641 
1642  // It is possible that an NT filename using backslashes is in the title string
1643  // this can give fits to unix codes where the backslash is assumed to be an escape
1644  // sequence. For the exodus file, just flip the backslash to a slash to prevent
1645  // this problem
1646 
1647  // Get a working copy of the title_string;
1648 
1649  char working_title[80];
1650  strncpy( working_title, title_string, 79 );
1651 
1652  int length = strlen( working_title );
1653  for( int pos = 0; pos < length; pos++ )
1654  {
1655  if( working_title[pos] == '\\' ) working_title[pos] = '/';
1656  }
1657 
1658  if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "title", length, working_title ) )
1659  {
1660  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define title attribute" );
1661  }
1662 
1663  // Add other attributes while we're at it
1664  float dum_vers = 6.28F;
1665  if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "api_version", NC_FLOAT, 1, &dum_vers ) )
1666  {
1667  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define api_version attribute" );
1668  }
1669  dum_vers = 6.28F;
1670  if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "version", NC_FLOAT, 1, &dum_vers ) )
1671  {
1672  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define version attribute" );
1673  }
1674  int dum_siz = sizeof( double );
1675  if( NC_NOERR != nc_put_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &dum_siz ) )
1676  {
1677  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define floating pt word size attribute" );
1678  }
1679 
1680  // Set up number of dimensions
1681 
1682  int num_el_blk, num_elem, num_nodes, num_dim, num_fa_blk, num_faces;
1683  if( nc_def_dim( ncFile, "num_dim", (size_t)mesh_info.num_dim, &num_dim ) != NC_NOERR )
1684  {
1685  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dimensions" );
1686  }
1687 
1688  if( nc_def_dim( ncFile, "num_nodes", mesh_info.num_nodes, &num_nodes ) != NC_NOERR )
1689  {
1690  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
1691  }
1692 
1693  int num_nod_per_fa; // it is needed for polyhedron only; need to compute it (connectivity of
1694  // faces of polyhedra)
1695  if( mesh_info.polyhedronFaces.size() > 0 )
1696  if( nc_def_dim( ncFile, "num_faces", (int)mesh_info.polyhedronFaces.size(), &num_faces ) != NC_NOERR )
1697  {
1698  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
1699  }
1700 
1701  if( nc_def_dim( ncFile, "num_elem", mesh_info.num_elements, &num_elem ) != NC_NOERR )
1702  {
1703  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements" );
1704  }
1705 
1706  if( nc_def_dim( ncFile, "num_el_blk", mesh_info.num_elementblocks, &num_el_blk ) != NC_NOERR )
1707  {
1708  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of element blocks" );
1709  }
1710 
1711  /* ...and some variables */
1712 
1713  /* Element block id status array */
1714  int idstat = -1;
1715  if( NC_NOERR != nc_def_var( ncFile, "eb_status", NC_LONG, 1, &num_el_blk, &idstat ) )
1716  {
1717  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block status array" );
1718  }
1719 
1720  /* Element block id array */
1721 
1722  int idarr = -1;
1723  if( NC_NOERR != nc_def_var( ncFile, "eb_prop1", NC_LONG, 1, &num_el_blk, &idarr ) )
1724  {
1725  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block id array" );
1726  }
1727 
1728  /* store property name as attribute of property array variable */
1729  if( NC_NOERR != nc_put_att_text( ncFile, idarr, "name", strlen( "ID" ), "ID" ) )
1730  {
1731  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element block property name ID" );
1732  }
1733 
1734  // count how many are polyhedron blocks
1735  int num_fa_blocks = 0; //, num_polyh_blocks = 0;
1736  for( unsigned int i = 0; i < block_data.size(); i++ )
1737  {
1738  MaterialSetData& block = block_data[i];
1739  if( EXOII_POLYHEDRON == block.element_type )
1740  {
1741  num_fa_blocks++;
1742  // num_polyh_blocks++;
1743  }
1744  }
1745  if( 0 == this->repeat_face_blocks && num_fa_blocks > 1 ) num_fa_blocks = 1;
1746  char wname[CHAR_STR_LEN];
1747 
1748  if( num_fa_blocks > 0 )
1749  {
1750  /* face block id status array */
1751  if( nc_def_dim( ncFile, "num_fa_blk", num_fa_blocks, &num_fa_blk ) != NC_NOERR )
1752  {
1753  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of face blocks" );
1754  }
1755 
1756  int idstatf = -1;
1757  if( NC_NOERR != nc_def_var( ncFile, "fa_status", NC_LONG, 1, &num_fa_blk, &idstatf ) )
1758  {
1759  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block status array" );
1760  }
1761 
1762  /* Element block id array */
1763 
1764  int idarrf = -1;
1765  if( NC_NOERR != nc_def_var( ncFile, "fa_prop1", NC_LONG, 1, &num_fa_blk, &idarrf ) )
1766  {
1767  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block id array" );
1768  }
1769 
1770  /* store property name as attribute of property array variable */
1771  if( NC_NOERR != nc_put_att_text( ncFile, idarrf, "name", strlen( "ID" ), "ID" ) )
1772  {
1773  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store face block property name ID" );
1774  }
1775  // determine the number of num_nod_per_face
1776  /*
1777  num_fa_in_blk1 = 15 ;
1778  num_nod_per_fa1 = 58 ;
1779 
1780  int fbconn1(num_nod_per_fa1) ;
1781  fbconn1:elem_type = "nsided" ;
1782  int fbepecnt1(num_fa_in_blk1) ;
1783  fbepecnt1:entity_type1 = "NODE" ;
1784  fbepecnt1:entity_type2 = "FACE" ;
1785  */
1786 
1787  int num_nodes_per_face = 0;
1788 
1789  int dims[1]; // maybe 1 is enough here
1790  for( Range::iterator eit = mesh_info.polyhedronFaces.begin(); eit != mesh_info.polyhedronFaces.end(); eit++ )
1791  {
1792  EntityHandle polyg = *eit;
1793  int nnodes = 0;
1794  const EntityHandle* conn = NULL;
1795  MB_CHK_ERR( mdbImpl->get_connectivity( polyg, conn, nnodes ) );
1796  num_nodes_per_face += nnodes;
1797  }
1798 
1799  // duplicate if needed; default is not duplicate
1800  for( int j = 1; j <= num_fa_blocks; j++ )
1801  {
1802  INS_ID( wname, "num_nod_per_fa%d", j );
1803  if( nc_def_dim( ncFile, wname, (size_t)num_nodes_per_face, &num_nod_per_fa ) != NC_NOERR )
1804  {
1805  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
1806  }
1807  dims[0] = num_nod_per_fa;
1808  INS_ID( wname, "fbconn%d", j ); // first one, or more
1809  int fbconn;
1810  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbconn ) )
1811  {
1812  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for face block " << 1 );
1813  }
1814  std::string element_type_string( "nsided" );
1815  if( NC_NOERR != nc_put_att_text( ncFile, fbconn, "elem_type", element_type_string.length(),
1816  element_type_string.c_str() ) )
1817  {
1818  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type nsided " );
1819  }
1820 
1821  INS_ID( wname, "num_fa_in_blk%d", j ); // first one, or more
1822  int num_fa_in_blk;
1823  if( nc_def_dim( ncFile, wname, (size_t)mesh_info.polyhedronFaces.size(), &num_fa_in_blk ) != NC_NOERR )
1824  {
1825  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
1826  }
1827 
1828  // fbepecnt
1829  INS_ID( wname, "fbepecnt%d", j ); // first one, or more
1830  int fbepecnt;
1831  dims[0] = num_fa_in_blk;
1832  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbepecnt ) )
1833  {
1834  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create fbepecnt array for block " << 1 );
1835  }
1836  std::string enttype1( "NODE" );
1837  if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type1", enttype1.length(), enttype1.c_str() ) )
1838  {
1839  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 1 " );
1840  }
1841  std::string enttype2( "FACE" );
1842  if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type2", enttype2.length(), enttype2.c_str() ) )
1843  {
1844  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 2 " );
1845  }
1846  }
1847  }
1848 
1849  // Define element blocks
1850 
1851  for( unsigned int i = 0; i < block_data.size(); i++ )
1852  {
1853  MaterialSetData& block = block_data[i];
1854 
1855  element_block_index = i + 1;
1856  int num_el_in_blk = -1, num_att_in_blk = -1;
1857  int blk_attrib, connect;
1858 
1859  /* Define number of elements in this block */
1860 
1861  INS_ID( wname, "num_el_in_blk%d", element_block_index );
1862  if( nc_def_dim( ncFile, wname, (size_t)block.number_elements, &num_el_in_blk ) != NC_NOERR )
1863  {
1864  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements/block for block " << i + 1 );
1865  }
1866 
1867  /* Define number of nodes per element for this block */
1868  INS_ID( wname, "num_nod_per_el%d", element_block_index );
1869  int num_nod_per_el = -1;
1870  if( EXOII_POLYHEDRON != block.element_type )
1871  if( nc_def_dim( ncFile, wname, (size_t)block.number_nodes_per_element, &num_nod_per_el ) != NC_NOERR )
1872  {
1873  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes/element for block " << block.id );
1874  }
1875 
1876  /* Define element attribute array for this block */
1877  int dims[3];
1878  if( block.number_attributes > 0 )
1879  {
1880  INS_ID( wname, "num_att_in_blk%d", element_block_index );
1881  if( nc_def_dim( ncFile, wname, (size_t)block.number_attributes, &num_att_in_blk ) != NC_NOERR )
1882  {
1883  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of attributes in block " << block.id );
1884  }
1885 
1886  INS_ID( wname, "attrib%d", element_block_index );
1887  dims[0] = num_el_in_blk;
1888  dims[1] = num_att_in_blk;
1889  if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 2, dims, &blk_attrib ) )
1890  {
1891  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define attributes for element block " << block.id );
1892  }
1893  }
1894 
1895  /* Define element connectivity array for this block */
1896 
1897  if( EXOII_POLYGON != block.element_type && EXOII_POLYHEDRON != block.element_type )
1898  {
1899  INS_ID( wname, "connect%d", element_block_index );
1900  dims[0] = num_el_in_blk;
1901  dims[1] = num_nod_per_el;
1902  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 2, dims, &connect ) )
1903  {
1904  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
1905  }
1906 
1907  /* Store element type as attribute of connectivity variable */
1908 
1909  std::string element_type_string( ExoIIUtil::ElementTypeNames[block.element_type] );
1910  if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1911  element_type_string.c_str() ) )
1912  {
1913  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
1914  }
1915  }
1916  else if( EXOII_POLYGON == block.element_type )
1917  {
1918  INS_ID( wname, "connect%d", element_block_index );
1919  // need to define num_nod_per_el as total number of nodes
1920  // ebepecnt1 as number of nodes per polygon
1921  /*
1922  * int connect1(num_nod_per_el1) ;
1923  connect1:elem_type = "nsided" ;
1924  int ebepecnt1(num_el_in_blk1) ;
1925  ebepecnt1:entity_type1 = "NODE" ;
1926  ebepecnt1:entity_type2 = "ELEM" ;*/
1927  dims[0] = num_nod_per_el;
1928  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &connect ) )
1929  {
1930  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
1931  }
1932  std::string element_type_string( "nsided" );
1933  if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
1934  element_type_string.c_str() ) )
1935  {
1936  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
1937  }
1938  INS_ID( wname, "ebepecnt%d", element_block_index );
1939  int ebepecnt;
1940  dims[0] = num_el_in_blk;
1941  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
1942  {
1943  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
1944  }
1945  std::string etype1( "NODE" );
1946  if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
1947  {
1948  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
1949  }
1950  std::string etype2( "ELEM" );
1951  if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
1952  {
1953  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
1954  }
1955  }
1956  else if( EXOII_POLYHEDRON == block.element_type )
1957  {
1958  // INS_ID(wname, "connect%d", element_block_index);
1959  /*
1960  testn face example: 3 polyh, 15 total faces, 2 shared; 15+2 = 17
1961  num_elem = 3 ;
1962  num_face = 15 ; // not needed?
1963  num_el_blk = 1 ;
1964 
1965  num_el_in_blk1 = 3 ;
1966  num_fac_per_el1 = 17 ;
1967 
1968  * num faces will be total face conn
1969  * num_faces_in_block will be number of faces (non repeated)
1970  * num_nodes_per_face will have total face connectivity
1971  */
1972  int num_faces2 = 0;
1973 
1974  for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
1975  {
1976  EntityHandle polyh = *eit;
1977  int nfaces = 0;
1978  const EntityHandle* conn = NULL;
1979  MB_CHK_ERR( mdbImpl->get_connectivity( polyh, conn, nfaces ) );
1980  num_faces2 += nfaces;
1981  }
1982 
1983  int num_fac_per_el;
1984 
1985  INS_ID( wname, "num_fac_per_el%d", element_block_index );
1986  if( nc_def_dim( ncFile, wname, (size_t)num_faces2, &num_fac_per_el ) != NC_NOERR )
1987  {
1988  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of faces per block " << block.id );
1989  }
1990 
1991  /*
1992  int facconn1(num_fac_per_el1) ;
1993  facconn1:elem_type = "NFACED" ;
1994  int ebepecnt1(num_el_in_blk1) ;
1995  ebepecnt1:entity_type1 = "FACE" ;
1996  ebepecnt1:entity_type2 = "ELEM" ;
1997  */
1998 
1999  // facconn
2000  INS_ID( wname, "facconn%d", element_block_index );
2001  int facconn;
2002  dims[0] = num_fac_per_el;
2003  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &facconn ) )
2004  {
2005  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create facconn array for block " << i + 1 );
2006  }
2007  std::string etype( "NFACED" );
2008  if( NC_NOERR != nc_put_att_text( ncFile, facconn, "elem_type", etype.length(), etype.c_str() ) )
2009  {
2010  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store elem type " << (int)block.element_type );
2011  }
2012 
2013  // ebepecnt
2014  INS_ID( wname, "ebepecnt%d", element_block_index );
2015  int ebepecnt;
2016  dims[0] = num_el_in_blk;
2017  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
2018  {
2019  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
2020  }
2021  std::string etype1( "FACE" );
2022  if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
2023  {
2024  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
2025  }
2026  std::string etype2( "ELEM" );
2027  if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
2028  {
2029  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
2030  }
2031 
2032  block.number_nodes_per_element = num_faces2; // connectivity for all polyhedra in block
2033  }
2034  }
2035 
2036  /* Node set id array: */
2037 
2038  int non_empty_nss = 0;
2039  // Need to go through nodesets to compute # nodesets, some might be empty
2040  std::vector< DirichletSetData >::iterator ns_it;
2041  for( ns_it = nodeset_data.begin(); ns_it != nodeset_data.end(); ++ns_it )
2042  {
2043  if( 0 != ( *ns_it ).number_nodes ) non_empty_nss++;
2044  }
2045 
2046  int num_ns = -1;
2047  int ns_idstat = -1, ns_idarr = -1;
2048  if( non_empty_nss > 0 )
2049  {
2050  if( nc_def_dim( ncFile, "num_node_sets", (size_t)( non_empty_nss ), &num_ns ) != NC_NOERR )
2051  {
2052  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of node sets" );
2053  }
2054 
2055  /* Node set id status array: */
2056 
2057  if( NC_NOERR != nc_def_var( ncFile, "ns_status", NC_LONG, 1, &num_ns, &ns_idstat ) )
2058  {
2059  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets status array" );
2060  }
2061 
2062  /* Node set id array: */
2063  if( NC_NOERR != nc_def_var( ncFile, "ns_prop1", NC_LONG, 1, &num_ns, &ns_idarr ) )
2064  {
2065  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets property array" );
2066  }
2067 
2068  /* Store property name as attribute of property array variable */
2069  if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "name", strlen( "ID" ), "ID" ) )
2070  {
2071  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store node set property name ID" );
2072  }
2073 
2074  // Now, define the arrays needed for each node set
2075 
2076  int index = 0;
2077 
2078  for( unsigned i = 0; i < nodeset_data.size(); i++ )
2079  {
2080  DirichletSetData node_set = nodeset_data[i];
2081 
2082  if( 0 == node_set.number_nodes )
2083  {
2084  MB_SET_ERR_CONT( "WriteNCDF: empty nodeset " << node_set.id );
2085  continue;
2086  }
2087  index++;
2088 
2089  int num_nod_ns = -1;
2090  INS_ID( wname, "num_nod_ns%d", index );
2091  if( nc_def_dim( ncFile, wname, (size_t)node_set.number_nodes, &num_nod_ns ) != NC_NOERR )
2092  {
2093  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for set " << node_set.id );
2094  }
2095 
2096  /* Create variable array in which to store the node set node list */
2097  int node_ns = -1;
2098  INS_ID( wname, "node_ns%d", index );
2099  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_nod_ns, &node_ns ) )
2100  {
2101  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node set " << node_set.id << " node list" );
2102  }
2103 
2104  // Create distribution factor array
2105  int fact_ns = -1;
2106  INS_ID( wname, "dist_fact_ns%d", index );
2107  if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 1, &num_nod_ns, &fact_ns ) )
2108  {
2109  MB_SET_ERR( MB_FAILURE,
2110  "WriteNCDF: failed to create node set " << node_set.id << " distribution factor list" );
2111  }
2112  }
2113  }
2114 
2115  /* Side set id array: */
2116 
2117  long non_empty_ss = 0;
2118  // Need to go through nodesets to compute # nodesets, some might be empty
2119  std::vector< NeumannSetData >::iterator ss_it;
2120  for( ss_it = sideset_data.begin(); ss_it != sideset_data.end(); ++ss_it )
2121  {
2122  if( 0 != ( *ss_it ).number_elements ) non_empty_ss++;
2123  }
2124 
2125  if( non_empty_ss > 0 )
2126  {
2127  int num_ss = -1;
2128  if( nc_def_dim( ncFile, "num_side_sets", non_empty_ss, &num_ss ) != NC_NOERR )
2129  {
2130  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of side sets" );
2131  }
2132 
2133  /* Side set id status array: */
2134  int ss_idstat = -1, ss_idarr = -1;
2135  if( NC_NOERR != nc_def_var( ncFile, "ss_status", NC_LONG, 1, &num_ss, &ss_idstat ) )
2136  {
2137  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set status" );
2138  }
2139 
2140  /* Side set id array: */
2141  if( NC_NOERR != nc_def_var( ncFile, "ss_prop1", NC_LONG, 1, &num_ss, &ss_idarr ) )
2142  {
2143  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set property" );
2144  }
2145 
2146  /* Store property name as attribute of property array variable */
2147  if( NC_NOERR != nc_put_att_text( ncFile, ss_idarr, "name", strlen( "ID" ), "ID" ) )
2148  {
2149  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store side set property name ID" );
2150  }
2151 
2152  // Now, define the arrays needed for each side set
2153 
2154  int index = 0;
2155  for( unsigned int i = 0; i < sideset_data.size(); i++ )
2156  {
2157  NeumannSetData side_set = sideset_data[i];
2158 
2159  // Don't define an empty set
2160  if( 0 == side_set.number_elements ) continue;
2161 
2162  index++;
2163 
2164  int num_side_ss = -1;
2165  int elem_ss = -1, side_ss = -1;
2166  INS_ID( wname, "num_side_ss%d", index );
2167  if( nc_def_dim( ncFile, wname, (size_t)side_set.number_elements, &num_side_ss ) != NC_NOERR )
2168  {
2169  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of sides in side set " << side_set.id );
2170  }
2171 
2172  INS_ID( wname, "elem_ss%d", index );
2173  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &elem_ss ) )
2174  {
2175  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element list for side set "
2176  << side_set.id ); /* Exit define mode and return */
2177  }
2178  INS_ID( wname, "side_ss%d", index );
2179  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &side_ss ) )
2180  {
2181  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create side list for side set "
2182  << side_set.id ); /* Exit define mode and return */
2183  }
2184 
2185  // sideset distribution factors
2186  int num_df_ss = -1;
2187  INS_ID( wname, "num_df_ss%d", index );
2188  if( nc_def_dim( ncFile, wname, (size_t)side_set.ss_dist_factors.size(), &num_df_ss ) != NC_NOERR )
2189  {
2190  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dist factors in side set "
2191  << side_set.id ); /* Exit define mode and return */
2192  }
2193 
2194  /* Create variable array in which to store the side set distribution factors */
2195 
2196  int fact_ss = -1;
2197  INS_ID( wname, "dist_fact_ss%d", index );
2198  if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_df_ss, &fact_ss ) )
2199  {
2200  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create dist factors list for side set "
2201  << side_set.id ); /* Exit define mode and return */
2202  }
2203  }
2204  }
2205 
2206  /* Node coordinate arrays: */
2207 
2208  int coord, name_coord, dims[3];
2209  dims[0] = num_dim;
2210  dims[1] = num_nodes;
2211  if( NC_NOERR != nc_def_var( ncFile, "coord", NC_DOUBLE, 2, dims, &coord ) )
2212  {
2213  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define node coordinate array" );
2214  }
2215 
2216  /* Coordinate names array */
2217 
2218  dims[0] = num_dim;
2219  dims[1] = dim_str;
2220  if( NC_NOERR != nc_def_var( ncFile, "coor_names", NC_CHAR, 2, dims, &name_coord ) )
2221  {
2222  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define coordinate name array" );
2223  }
2224 
2225  // Define genesis maps if required
2226 
2227  if( write_maps )
2228  {
2229  // Element map
2230  int elem_map = -1, elem_map2 = -1, node_map = -1;
2231  if( NC_NOERR != nc_def_var( ncFile, "elem_map", NC_LONG, 1, &num_elem, &elem_map ) )
2232  {
2233  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element map array" ); /* Exit define mode and return */
2234  }
2235 
2236  // Create the element number map
2237  if( NC_NOERR != nc_def_var( ncFile, "elem_num_map", NC_LONG, 1, &num_elem, &elem_map2 ) )
2238  {
2239  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element numbering map" ); /* Exit define mode
2240  and return */
2241  }
2242 
2243  // Create node number map
2244  if( NC_NOERR != nc_def_var( ncFile, "node_num_map", NC_LONG, 1, &num_nodes, &node_map ) )
2245  {
2246  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node numbering map array" ); /* Exit define mode and
2247  return */
2248  }
2249  }
2250 
2251  // Define qa records to be used
2252 
2253  int num_qa_rec = mesh_info.qaRecords.size() / 4;
2254  int num_qa = -1;
2255 
2256  if( nc_def_dim( ncFile, "num_qa_rec", (long)num_qa_rec, &num_qa ) != NC_NOERR )
2257  {
2258  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array size" );
2259  }
2260 
2261  // Define qa array
2262  int qa_title;
2263  dims[0] = num_qa;
2264  dims[1] = dim_four;
2265  dims[2] = dim_str;
2266  if( NC_NOERR != nc_def_var( ncFile, "qa_records", NC_CHAR, 3, dims, &qa_title ) )
2267  {
2268  MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array" );
2269  }
2270 
2271  // Take it out of define mode
2272  if( NC_NOERR != nc_enddef( ncFile ) )
2273  {
2274  MB_SET_ERR( MB_FAILURE, "WriteNCDF: Trouble leaving define mode" );
2275  }
2276 
2277  return MB_SUCCESS;
2278 }
2279 
2280 ErrorCode WriteNCDF::open_file( const char* filename )
2281 {
2282  // Not a valid filename
2283  if( strlen( (const char*)filename ) == 0 )
2284  {
2285  MB_SET_ERR( MB_FAILURE, "Output Exodus filename not specified" );
2286  }
2287 
2288  int fail = nc_create( filename, NC_CLOBBER, &ncFile );
2289 
2290  // File couldn't be opened
2291  if( NC_NOERR != fail )
2292  {
2293  MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
2294  }
2295 
2296  return MB_SUCCESS;
2297 }
2298 
2299 ErrorCode WriteNCDF::get_sideset_elems( EntityHandle sideset,
2300  int current_sense,
2301  Range& forward_elems,
2302  Range& reverse_elems )
2303 {
2304  Range ss_elems, ss_meshsets;
2305 
2306  // Get the sense tag; don't need to check return, might be an error if the tag
2307  // hasn't been created yet
2308  Tag sense_tag = 0;
2309  mdbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
2310 
2311  // Get the entities in this set
2312  ErrorCode result = mdbImpl->get_entities_by_handle( sideset, ss_elems, true );
2313  if( MB_FAILURE == result ) return result;
2314 
2315  // Now remove the meshsets into the ss_meshsets; first find the first meshset,
2316  Range::iterator range_iter = ss_elems.begin();
2317  while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
2318  ++range_iter;
2319 
2320  // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
2321  if( range_iter != ss_elems.end() )
2322  {
2323  std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
2324  ss_elems.erase( range_iter, ss_elems.end() );
2325  }
2326 
2327  // OK, for the elements, check the sense of this set and copy into the right range
2328  // (if the sense is 0, copy into both ranges)
2329 
2330  // Need to step forward on list until we reach the right dimension
2331  Range::iterator dum_it = ss_elems.end();
2332  --dum_it;
2333  int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
2334  dum_it = ss_elems.begin();
2335  while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
2336  ++dum_it;
2337 
2338  if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
2339  if( current_sense == -1 || current_sense == 0 )
2340  std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
2341 
2342  // Now loop over the contained meshsets, getting the sense of those and calling this
2343  // function recursively
2344  for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
2345  {
2346  // First get the sense; if it's not there, by convention it's forward
2347  int this_sense;
2348  if( 0 == sense_tag || MB_FAILURE == mdbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
2349  this_sense = 1;
2350 
2351  // Now get all the entities on this meshset, with the proper (possibly reversed) sense
2352  get_sideset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
2353  }
2354 
2355  return result;
2356 }
2357 
2358 } // namespace moab