Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
WriteCCMIO.cpp
Go to the documentation of this file.
1 /*
2  * CCMIO file structure
3  *
4  * Root
5  * State(kCCMIOState)
6  * Processor*
7  * VerticesID
8  * TopologyID
9  * InitialID
10  * SolutionID
11  * Vertices*
12  * ->WriteVerticesx, WriteMap
13  * Topology*
14  * Boundary faces*(kCCMIOBoundaryFaces)
15  * ->WriteFaces, WriteFaceCells, WriteMap
16  * Internal faces(kCCMIOInternalFaces)
17  * Cells (kCCMIOCells)
18  * ->WriteCells (mapID), WriteMap, WriteCells
19  * Solution
20  * Phase
21  * Field
22  * FieldData
23  * Problem(kCCMIOProblemDescription)
24  * CellType* (kCCMIOCellType)
25  * Index (GetEntityIndex), MaterialId(WriteOpti), MaterialType(WriteOptstr),
26  * PorosityId(WriteOpti), SpinId(WriteOpti), GroupId(WriteOpti)
27  *
28  * MaterialType (CCMIOWriteOptstr in readexample)
29  * constants (see readexample)
30  * lagrangian data (CCMIOWriteLagrangianData)
31  * vertices label (CCMIOEntityDescription)
32  * restart info: char solver[], iterations, time, char timeUnits[], angle
33  * (CCMIOWriteRestartInfo, kCCMIORestartData), reference data?
34  * phase:
35  * field: char name[], dims, CCMIODataType datatype, char units[]
36  * dims = kCCMIOScalar (CCMIOWriteFieldDataf),
37  * kCCMIOVector (CCMIOWriteMultiDimensionalFieldData),
38  * kCCMIOTensor
39  * MonitoringSets: num, name (CellSet, VertexSet, BoundarySet, BlockSet, SplineSet, CoupleSet)
40  * CCMIOGetProstarSet, CCMIOWriteOpt1i,
41  */
42 
43 #ifdef WIN32
44 #ifdef _DEBUG
45 // turn off warnings that say they debugging identifier has been truncated
46 // this warning comes up when using some STL containers
47 #pragma warning( disable : 4786 )
48 #endif
49 #endif
50 
51 #include "WriteCCMIO.hpp"
52 #include "ccmio.h"
53 #include "ccmioutility.h"
54 #include "ccmiocore.h"
55 #include <utility>
56 #include <algorithm>
57 #include <ctime>
58 #include <string>
59 #include <vector>
60 #include <cstdio>
61 #include <iostream>
62 #include <algorithm>
63 #include <sstream>
64 
65 #include "moab/Interface.hpp"
66 #include "moab/Range.hpp"
67 #include "moab/CN.hpp"
68 #include "moab/Skinner.hpp"
69 #include <cassert>
70 #include "Internals.hpp"
71 #include "ExoIIUtil.hpp"
72 #include "MBTagConventions.hpp"
73 #ifdef MOAB_HAVE_MPI
74 #include "MBParallelConventions.h"
75 #endif
76 #include "moab/WriteUtilIface.hpp"
77 
78 namespace moab
79 {
80 
81 static char const kStateName[] = "default";
82 
83 /*
84  static const int ccm_types[] = {
85  1, // MBVERTEX
86  2, // MBEDGE
87  -1, // MBTRI
88  -1, // MBQUAD
89  -1, // MBPOLYGON
90  13, // MBTET
91  14, // MBPYRAMID
92  12, // MBPRISM
93  -1, // MBKNIFE
94  11, // MBHEX
95  255 // MBPOLYHEDRON
96  };
97 */
98 
99 #define CHK_SET_CCMERR( ccm_err_code, ccm_err_msg ) \
100  { \
101  if( kCCMIONoErr != ( ccm_err_code ) ) MB_SET_ERR( MB_FAILURE, ccm_err_msg ); \
102  }
103 
105 {
106  return new WriteCCMIO( iface );
107 }
108 
110  : mbImpl( impl ), mCurrentMeshHandle( 0 ), mPartitionSetTag( 0 ), mNameTag( 0 ), mMaterialIdTag( 0 ),
111  mMaterialTypeTag( 0 ), mRadiationTag( 0 ), mPorosityIdTag( 0 ), mSpinIdTag( 0 ), mGroupIdTag( 0 ),
112  mColorIdxTag( 0 ), mProcessorIdTag( 0 ), mLightMaterialTag( 0 ), mFreeSurfaceMaterialTag( 0 ), mThicknessTag( 0 ),
113  mProstarRegionNumberTag( 0 ), mBoundaryTypeTag( 0 ), mCreatingProgramTag( 0 ), mDimension( 0 ),
114  mWholeMesh( false )
115 {
116  assert( impl != NULL );
117 
118  impl->query_interface( mWriteIface );
119 
120  // Initialize in case tag_get_handle fails below
121  //! Get and cache predefined tag handles
122  int negone = -1;
124  &negone );
125 
127  &negone );
128 
130  &negone );
131 
132  mGlobalIdTag = impl->globalId_tag();
133 
134 #ifdef MOAB_HAVE_MPI
136  // No need to check result, if it's not there, we don't create one
137 #endif
138 
139  int dum_val_array[] = { -1, -1, -1, -1 };
141  dum_val_array );
142 
143  impl->tag_get_handle( "__WriteCCMIO element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
144 
145  // Don't need to check return of following, since it doesn't matter if there isn't one
147 }
148 
150 {
153 }
154 
155 ErrorCode WriteCCMIO::write_file( const char* file_name,
156  const bool overwrite,
157  const FileOptions&,
158  const EntityHandle* ent_handles,
159  const int num_sets,
160  const std::vector< std::string >& /* qa_list */,
161  const Tag* /* tag_list */,
162  int /* num_tags */,
163  int /* export_dimension */ )
164 {
165  assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
166 
167  ErrorCode result;
168 
169  // Check overwrite flag and file existence
170  if( !overwrite )
171  {
172  FILE* file = fopen( file_name, "r" );
173  if( file )
174  {
175  fclose( file );
176  MB_SET_ERR( MB_FILE_WRITE_ERROR, "File exists but overwrite set to false" );
177  }
178  }
179 
180  mDimension = 3;
181 
182  std::vector< EntityHandle > matsets, dirsets, neusets, partsets;
183 
184  // Separate into material, dirichlet, neumann, partition sets
185  result = get_sets( ent_handles, num_sets, matsets, dirsets, neusets, partsets );MB_CHK_SET_ERR( result, "Failed to get material/etc. sets" );
186 
187  // If entity handles were input but didn't contain matsets, return error
188  if( ent_handles && matsets.empty() )
189  {
190  MB_SET_ERR( MB_FILE_WRITE_ERROR, "Sets input to write but no material sets found" );
191  }
192 
193  // Otherwise, if no matsets, use root set
194  if( matsets.empty() ) matsets.push_back( 0 );
195 
196  std::vector< MaterialSetData > matset_info;
197  Range all_verts;
198  result = gather_matset_info( matsets, matset_info, all_verts );MB_CHK_SET_ERR( result, "gathering matset info failed" );
199 
200  // Assign vertex gids
201  result = mWriteIface->assign_ids( all_verts, mGlobalIdTag, 1 );MB_CHK_SET_ERR( result, "Failed to assign vertex global ids" );
202 
203  // Some CCMIO descriptors
204  CCMIOID rootID, topologyID, stateID, problemID, verticesID, processorID;
205 
206  // Try to open the file and establish state
207  result = open_file( file_name, overwrite, rootID );MB_CHK_SET_ERR( result, "Couldn't open file or create state" );
208 
209  result = create_ccmio_structure( rootID, stateID, processorID );MB_CHK_SET_ERR( result, "Problem creating CCMIO file structure" );
210 
211  result = write_nodes( rootID, all_verts, mDimension, verticesID );MB_CHK_SET_ERR( result, "write_nodes failed" );
212 
213  std::vector< NeumannSetData > neuset_info;
214  result = gather_neuset_info( neusets, neuset_info );MB_CHK_SET_ERR( result, "Failed to get neumann set info" );
215 
216  result = write_cells_and_faces( rootID, matset_info, neuset_info, all_verts, topologyID );MB_CHK_SET_ERR( result, "write_cells_and_faces failed" );
217 
218  result = write_problem_description( rootID, stateID, problemID, processorID, matset_info, neuset_info );MB_CHK_SET_ERR( result, "write_problem_description failed" );
219 
220  result = write_solution_data();MB_CHK_SET_ERR( result, "Trouble writing solution data" );
221 
222  result = write_processor( processorID, verticesID, topologyID );MB_CHK_SET_ERR( result, "Trouble writing processor" );
223 
224  result = close_and_compress( file_name, rootID );MB_CHK_SET_ERR( result, "Close or compress failed" );
225 
226  return MB_SUCCESS;
227 }
228 
230 {
231  // For now, no solution (tag) data
232  return MB_SUCCESS;
233 }
234 
235 ErrorCode WriteCCMIO::write_processor( CCMIOID processorID, CCMIOID verticesID, CCMIOID topologyID )
236 {
237  CCMIOError error = kCCMIONoErr;
238 
239  // Now we have the mesh (vertices and topology) and the post data written.
240  // Since we now have their IDs, we can write out the processor information.
241  CCMIOWriteProcessor( &error, processorID, NULL, &verticesID, NULL, &topologyID, NULL, NULL, NULL, NULL );
242  CHK_SET_CCMERR( error, "Problem writing CCMIO processor" );
243 
244  return MB_SUCCESS;
245 }
246 
247 ErrorCode WriteCCMIO::create_ccmio_structure( CCMIOID rootID, CCMIOID& stateID, CCMIOID& processorID )
248 {
249  // Create problem state and other CCMIO nodes under it
250  CCMIOError error = kCCMIONoErr;
251 
252  // Create a new state (or re-use an existing one).
253  if( CCMIOGetState( NULL, rootID, kStateName, NULL, &stateID ) != kCCMIONoErr )
254  {
255  CCMIONewState( &error, rootID, kStateName, NULL, NULL, &stateID );
256  CHK_SET_CCMERR( error, "Trouble creating state" );
257  }
258 
259  // Create or get an old processor for this state
260  CCMIOSize_t i = CCMIOSIZEC( 0 );
261  if( CCMIONextEntity( NULL, stateID, kCCMIOProcessor, &i, &processorID ) != kCCMIONoErr )
262  {
263  CCMIONewEntity( &error, stateID, kCCMIOProcessor, NULL, &processorID );
264  CHK_SET_CCMERR( error, "Trouble creating processor node" );
265  }
266  // Get rid of any data that may be in this processor (if the state was
267  // not new).
268  else
269  {
270  CCMIOClearProcessor( &error, stateID, processorID, TRUE, TRUE, TRUE, TRUE, TRUE );
271  CHK_SET_CCMERR( error, "Trouble clearing processor data" );
272  }
273 
274  /*
275  // for (; i < CCMIOSIZEC(partsets.size()); i++) {
276  CCMIOSize_t id = CCMIOSIZEC(0);
277  if (CCMIONextEntity(NULL, stateID, kCCMIOProcessor, &id, &processorID) != kCCMIONoErr)
278  CCMIONewEntity(&error, stateID, kCCMIOProcessor, NULL, &processorID);
279  CHKCCMERR(error, "Trouble creating processor node.");
280  */
281  return MB_SUCCESS;
282 }
283 
284 ErrorCode WriteCCMIO::close_and_compress( const char*, CCMIOID rootID )
285 {
286  CCMIOError error = kCCMIONoErr;
287  CCMIOCloseFile( &error, rootID );
288  CHK_SET_CCMERR( error, "File close failed" );
289 
290  // The CCMIO library uses ADF to store the actual data. Unfortunately,
291  // ADF leaks disk space; deleting a node does not recover all the disk
292  // space. Now that everything is successfully written it might be useful
293  // to call CCMIOCompress() here to ensure that the file is as small as
294  // possible. Please see the Core API documentation for caveats on its
295  // usage.
296  // CCMIOCompress(&error, const_cast<char*>(filename));CHK_SET_CCMERR(error, "Error compressing
297  // file");
298 
299  return MB_SUCCESS;
300 }
301 
302 ErrorCode WriteCCMIO::open_file( const char* filename, bool, CCMIOID& rootID )
303 {
304  CCMIOError error = kCCMIONoErr;
305  CCMIOOpenFile( &error, filename, kCCMIOWrite, &rootID );
306  CHK_SET_CCMERR( error, "Cannot open file" );
307 
308  return MB_SUCCESS;
309 }
310 
312  int num_sets,
313  std::vector< EntityHandle >& matsets,
314  std::vector< EntityHandle >& dirsets,
315  std::vector< EntityHandle >& neusets,
316  std::vector< EntityHandle >& partsets )
317 {
318  if( num_sets == 0 )
319  {
320  // Default to all defined sets
321  Range this_range;
322  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
323  std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
324  this_range.clear();
325  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
326  std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
327  this_range.clear();
328  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
329  std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
330  if( mPartitionSetTag )
331  {
332  this_range.clear();
333  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mPartitionSetTag, NULL, 1, this_range );
334  std::copy( this_range.begin(), this_range.end(), std::back_inserter( partsets ) );
335  }
336  }
337  else
338  {
339  int dummy;
340  for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
341  {
342  if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
343  matsets.push_back( *iter );
344  else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
345  dirsets.push_back( *iter );
346  else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
347  neusets.push_back( *iter );
348  else if( mPartitionSetTag && MB_SUCCESS == mbImpl->tag_get_data( mPartitionSetTag, &( *iter ), 1, &dummy ) )
349  partsets.push_back( *iter );
350  }
351  }
352 
353  return MB_SUCCESS;
354 }
355 
357  CCMIOID stateID,
358  CCMIOID& problemID,
359  CCMIOID processorID,
360  std::vector< WriteCCMIO::MaterialSetData >& matset_data,
361  std::vector< WriteCCMIO::NeumannSetData >& neuset_data )
362 {
363  // Write out a dummy problem description. If we happen to know that
364  // there already is a problem description previously recorded that
365  // is valid we could skip this step.
366  CCMIOID id;
367  CCMIOError error = kCCMIONoErr;
368  ErrorCode rval;
369  const EntityHandle mesh = 0;
370 
371  bool root_tagged = false, other_set_tagged = false;
372  Tag simname;
373  Range dum_sets;
374  rval = mbImpl->tag_get_handle( "Title", 0, MB_TYPE_OPAQUE, simname, MB_TAG_ANY );
375  if( MB_SUCCESS == rval )
376  {
377  int tag_size;
378  rval = mbImpl->tag_get_bytes( simname, tag_size );
379  if( MB_SUCCESS == rval )
380  {
381  std::vector< char > title_tag( tag_size + 1 );
382  rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &simname, NULL, 1, dum_sets );
383  if( MB_SUCCESS == rval && !dum_sets.empty() )
384  {
385  rval = mbImpl->tag_get_data( simname, &( *dum_sets.begin() ), 1, &title_tag[0] );MB_CHK_SET_ERR( rval, "Problem getting simulation name tag" );
386  other_set_tagged = true;
387  }
388  else if( MB_SUCCESS == rval )
389  {
390  // Check to see if interface was tagged
391  rval = mbImpl->tag_get_data( simname, &mesh, 1, &title_tag[0] );
392  if( MB_SUCCESS == rval )
393  root_tagged = true;
394  else
395  rval = MB_SUCCESS;
396  }
397  *title_tag.rbegin() = '\0';
398  if( root_tagged || other_set_tagged )
399  {
400  CCMIONode rootNode;
401  if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
402  {
403  CCMIOSetTitle( &error, rootNode, &title_tag[0] );
404  CHK_SET_CCMERR( error, "Trouble setting title" );
405  }
406  }
407  }
408  }
409 
410  rval = mbImpl->tag_get_handle( "CreatingProgram", 0, MB_TYPE_OPAQUE, mCreatingProgramTag, MB_TAG_ANY );
411  if( MB_SUCCESS == rval )
412  {
413  int tag_size;
414  rval = mbImpl->tag_get_bytes( mCreatingProgramTag, tag_size );
415  if( MB_SUCCESS == rval )
416  {
417  std::vector< char > cp_tag( tag_size + 1 );
418  rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mCreatingProgramTag, NULL, 1, dum_sets );
419  if( MB_SUCCESS == rval && !dum_sets.empty() )
420  {
421  rval = mbImpl->tag_get_data( mCreatingProgramTag, &( *dum_sets.begin() ), 1, &cp_tag[0] );MB_CHK_SET_ERR( rval, "Problem getting creating program tag" );
422  other_set_tagged = true;
423  }
424  else if( MB_SUCCESS == rval )
425  {
426  // Check to see if interface was tagged
427  rval = mbImpl->tag_get_data( mCreatingProgramTag, &mesh, 1, &cp_tag[0] );
428  if( MB_SUCCESS == rval )
429  root_tagged = true;
430  else
431  rval = MB_SUCCESS;
432  }
433  *cp_tag.rbegin() = '\0';
434  if( root_tagged || other_set_tagged )
435  {
436  CCMIONode rootNode;
437  if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
438  {
439  CCMIOWriteOptstr( &error, processorID, "CreatingProgram", &cp_tag[0] );
440  CHK_SET_CCMERR( error, "Trouble setting creating program" );
441  }
442  }
443  }
444  }
445 
446  CCMIONewEntity( &error, rootID, kCCMIOProblemDescription, NULL, &problemID );
447  CHK_SET_CCMERR( error, "Trouble creating problem node" );
448 
449  // Write material types and other info
450  for( unsigned int i = 0; i < matset_data.size(); i++ )
451  {
452  if( !matset_data[i].setName.empty() )
453  {
454  CCMIONewIndexedEntity( &error, problemID, kCCMIOCellType, matset_data[i].matsetId,
455  matset_data[i].setName.c_str(), &id );
456  CHK_SET_CCMERR( error, "Failure creating celltype node" );
457 
458  CCMIOWriteOptstr( &error, id, "MaterialType", matset_data[i].setName.c_str() );
459  CHK_SET_CCMERR( error, "Error assigning material name" );
460  }
461  else
462  {
463  char dum_name[NAME_TAG_SIZE];
464  std::ostringstream os;
465  std::string mat_name = "Material", temp_str;
466  os << mat_name << ( i + 1 );
467  temp_str = os.str();
468  strcpy( dum_name, temp_str.c_str() );
469  CCMIONewIndexedEntity( &error, problemID, kCCMIOCellType, matset_data[i].matsetId, dum_name, &id );
470  CHK_SET_CCMERR( error, "Failure creating celltype node" );
471 
472  CCMIOWriteOptstr( &error, id, "MaterialType", dum_name );
473  CHK_SET_CCMERR( error, "Error assigning material name" );
474 
475  os.str( "" );
476  }
477  rval = write_int_option( "MaterialId", matset_data[i].setHandle, mMaterialIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing MaterialId option" );
478 
479  rval = write_int_option( "Radiation", matset_data[i].setHandle, mRadiationTag, id );MB_CHK_SET_ERR( rval, "Trouble writing Radiation option" );
480 
481  rval = write_int_option( "PorosityId", matset_data[i].setHandle, mPorosityIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing PorosityId option" );
482 
483  rval = write_int_option( "SpinId", matset_data[i].setHandle, mSpinIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing SpinId option" );
484 
485  rval = write_int_option( "GroupId", matset_data[i].setHandle, mGroupIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing GroupId option" );
486 
487  rval = write_int_option( "ColorIdx", matset_data[i].setHandle, mColorIdxTag, id );MB_CHK_SET_ERR( rval, "Trouble writing ColorIdx option" );
488 
489  rval = write_int_option( "ProcessorId", matset_data[i].setHandle, mProcessorIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing ProcessorId option" );
490 
491  rval = write_int_option( "LightMaterial", matset_data[i].setHandle, mLightMaterialTag, id );MB_CHK_SET_ERR( rval, "Trouble writing LightMaterial option." );
492 
493  rval = write_int_option( "FreeSurfaceMaterial", matset_data[i].setHandle, mFreeSurfaceMaterialTag, id );MB_CHK_SET_ERR( rval, "Trouble writing FreeSurfaceMaterial option" );
494 
495  rval = write_dbl_option( "Thickness", matset_data[i].setHandle, mThicknessTag, id );MB_CHK_SET_ERR( rval, "Trouble writing Thickness option" );
496 
497  rval = write_str_option( "MaterialType", matset_data[i].setHandle, mMaterialTypeTag, id );MB_CHK_SET_ERR( rval, "Trouble writing MaterialType option" );
498  }
499 
500  // Write neumann set info
501  for( unsigned int i = 0; i < neuset_data.size(); i++ )
502  {
503  // Use the label to encode the id
504  std::ostringstream dum_id;
505  dum_id << neuset_data[i].neusetId;
506  CCMIONewIndexedEntity( &error, problemID, kCCMIOBoundaryRegion, neuset_data[i].neusetId, dum_id.str().c_str(),
507  &id );
508  CHK_SET_CCMERR( error, "Failure creating BoundaryRegion node" );
509 
510  rval = write_str_option( "BoundaryName", neuset_data[i].setHandle, mNameTag, id );MB_CHK_SET_ERR( rval, "Trouble writing boundary type number" );
511 
512  rval = write_str_option( "BoundaryType", neuset_data[i].setHandle, mBoundaryTypeTag, id );MB_CHK_SET_ERR( rval, "Trouble writing boundary type number" );
513 
514  rval = write_int_option( "ProstarRegionNumber", neuset_data[i].setHandle, mProstarRegionNumberTag, id );MB_CHK_SET_ERR( rval, "Trouble writing prostar region number" );
515  }
516 
517  CCMIOWriteState( &error, stateID, problemID, "Example state" );
518  CHK_SET_CCMERR( error, "Failure writing problem state" );
519 
520  // Get cell types; reuse cell ids array
521  // for (i = 0, rit = all_elems.begin(); i < num_elems; i++, ++rit) {
522  // egids[i] = ccm_types[mbImpl->type_from_handle(*rit)];
523  // assert(-1 != egids[i]);
524  // }
525 
526  return MB_SUCCESS;
527 }
528 
529 ErrorCode WriteCCMIO::write_int_option( const char* opt_name, EntityHandle seth, Tag& tag, CCMIOID& node )
530 {
531  ErrorCode rval;
532 
533  if( !tag )
534  {
535  rval = mbImpl->tag_get_handle( opt_name, 1, MB_TYPE_INTEGER, tag );
536  // Return success since that just means we don't have to write this option
537  if( MB_SUCCESS != rval ) return MB_SUCCESS;
538  }
539 
540  int dum_val;
541  rval = mbImpl->tag_get_data( tag, &seth, 1, &dum_val );
542  // Return success since that just means we don't have to write this option
543  if( MB_SUCCESS != rval ) return MB_SUCCESS;
544 
545  CCMIOError error = kCCMIONoErr;
546  CCMIOWriteOpti( &error, node, opt_name, dum_val );
547  CHK_SET_CCMERR( error, "Trouble writing int option" );
548 
549  return MB_SUCCESS;
550 }
551 
552 ErrorCode WriteCCMIO::write_dbl_option( const char* opt_name, EntityHandle seth, Tag& tag, CCMIOID& node )
553 {
554  ErrorCode rval;
555 
556  if( !tag )
557  {
558  rval = mbImpl->tag_get_handle( opt_name, 1, MB_TYPE_DOUBLE, tag );
559  // Return success since that just means we don't have to write this option
560  if( MB_SUCCESS != rval ) return MB_SUCCESS;
561  }
562 
563  double dum_val;
564  rval = mbImpl->tag_get_data( tag, &seth, 1, &dum_val );
565  // Return success since that just means we don't have to write this option
566  if( MB_SUCCESS != rval ) return MB_SUCCESS;
567 
568  CCMIOError error = kCCMIONoErr;
569  CCMIOWriteOptf( &error, node, opt_name, dum_val );
570  CHK_SET_CCMERR( error, "Trouble writing int option" );
571 
572  return MB_SUCCESS;
573 }
574 
576  EntityHandle seth,
577  Tag& tag,
578  CCMIOID& node,
579  const char* other_name )
580 {
581  int tag_size;
582  ErrorCode rval;
583 
584  if( !tag )
585  {
586  rval = mbImpl->tag_get_handle( opt_name, 0, MB_TYPE_OPAQUE, tag, MB_TAG_ANY );
587  // Return success since that just means we don't have to write this option
588  if( MB_SUCCESS != rval ) return MB_SUCCESS;
589  }
590 
591  rval = mbImpl->tag_get_bytes( tag, tag_size );
592  if( MB_SUCCESS != rval ) return MB_SUCCESS;
593  std::vector< char > opt_val( tag_size + 1 );
594 
595  rval = mbImpl->tag_get_data( tag, &seth, 1, &opt_val[0] );
596  if( MB_SUCCESS != rval ) return MB_SUCCESS;
597 
598  // Null-terminate if necessary
599  if( std::find( opt_val.begin(), opt_val.end(), '\0' ) == opt_val.end() ) *opt_val.rbegin() = '\0';
600 
601  CCMIOError error = kCCMIONoErr;
602  if( other_name )
603  {
604  CCMIOWriteOptstr( &error, node, other_name, &opt_val[0] );
605  CHK_SET_CCMERR( error, "Failure writing an option string MaterialType" );
606  }
607  else
608  {
609  CCMIOWriteOptstr( &error, node, opt_name, &opt_val[0] );
610  CHK_SET_CCMERR( error, "Failure writing an option string MaterialType" );
611  }
612 
613  return MB_SUCCESS;
614 }
615 
616 ErrorCode WriteCCMIO::gather_matset_info( std::vector< EntityHandle >& matsets,
617  std::vector< MaterialSetData >& matset_data,
618  Range& all_verts )
619 {
620  ErrorCode result;
621  matset_data.resize( matsets.size() );
622  if( 1 == matsets.size() && 0 == matsets[0] )
623  {
624  // Whole mesh
625  mWholeMesh = true;
626 
627  result = mbImpl->get_entities_by_dimension( 0, mDimension, matset_data[0].elems );MB_CHK_SET_ERR( result, "Trouble getting all elements in mesh" );
628  result = mWriteIface->gather_nodes_from_elements( matset_data[0].elems, mEntityMark, all_verts );MB_CHK_SET_ERR( result, "Trouble gathering nodes from elements" );
629 
630  return result;
631  }
632 
633  std::vector< unsigned char > marks;
634  for( unsigned int i = 0; i < matsets.size(); i++ )
635  {
636  EntityHandle this_set = matset_data[i].setHandle = matsets[i];
637 
638  // Get all Entity Handles in the set
639  result = mbImpl->get_entities_by_dimension( this_set, mDimension, matset_data[i].elems, true );MB_CHK_SET_ERR( result, "Trouble getting m-dimensional ents" );
640 
641  // Get all connected vertices
642  result = mWriteIface->gather_nodes_from_elements( matset_data[i].elems, mEntityMark, all_verts );MB_CHK_SET_ERR( result, "Trouble getting vertices for a matset" );
643 
644  // Check for consistent entity type
645  EntityType start_type = mbImpl->type_from_handle( *matset_data[i].elems.begin() );
646  if( start_type == mbImpl->type_from_handle( *matset_data[i].elems.rbegin() ) )
647  matset_data[i].entityType = start_type;
648 
649  // Mark elements in this matset
650  marks.resize( matset_data[i].elems.size(), 0x1 );
651  result = mbImpl->tag_set_data( mEntityMark, matset_data[i].elems, &marks[0] );MB_CHK_SET_ERR( result, "Couln't mark entities being output" );
652 
653  // Get id for this matset
654  result = mbImpl->tag_get_data( mMaterialSetTag, &this_set, 1, &matset_data[i].matsetId );MB_CHK_SET_ERR( result, "Couln't get global id for material set" );
655 
656  // Get name for this matset
657  if( mNameTag )
658  {
659  char dum_name[NAME_TAG_SIZE];
660  result = mbImpl->tag_get_data( mNameTag, &this_set, 1, dum_name );
661  if( MB_SUCCESS == result ) matset_data[i].setName = dum_name;
662 
663  // Reset success, so later checks don't fail
664  result = MB_SUCCESS;
665  }
666  }
667 
668  if( all_verts.empty() )
669  {
670  MB_SET_ERR( MB_FILE_WRITE_ERROR, "No vertices from elements" );
671  }
672 
673  return MB_SUCCESS;
674 }
675 
676 ErrorCode WriteCCMIO::gather_neuset_info( std::vector< EntityHandle >& neusets,
677  std::vector< NeumannSetData >& neuset_info )
678 {
679  ErrorCode result;
680 
681  neuset_info.resize( neusets.size() );
682  for( unsigned int i = 0; i < neusets.size(); i++ )
683  {
684  EntityHandle this_set = neuset_info[i].setHandle = neusets[i];
685 
686  // Get all Entity Handles of one less dimension than that being output
687  result = mbImpl->get_entities_by_dimension( this_set, mDimension - 1, neuset_info[i].elems, true );MB_CHK_SET_ERR( result, "Trouble getting (m-1)-dimensional ents for neuset" );
688 
689  result = mbImpl->tag_get_data( mGlobalIdTag, &this_set, 1, &neuset_info[i].neusetId );
690  if( MB_TAG_NOT_FOUND == result )
691  {
692  result = mbImpl->tag_get_data( mNeumannSetTag, &this_set, 1, &neuset_info[i].neusetId );
693  if( MB_SUCCESS != result )
694  // Need some id; use the loop iteration number
695  neuset_info[i].neusetId = i;
696  }
697 
698  // Get name for this neuset
699  if( mNameTag )
700  {
701  char dum_name[NAME_TAG_SIZE];
702  result = mbImpl->tag_get_data( mNameTag, &this_set, 1, dum_name );
703  if( MB_SUCCESS == result ) neuset_info[i].setName = dum_name;
704 
705  // Reset success, so later checks don't fail
706  result = MB_SUCCESS;
707  }
708  }
709 
710  return MB_SUCCESS;
711 }
712 
713 ErrorCode WriteCCMIO::get_gids( const Range& ents, int*& gids, int& minid, int& maxid )
714 {
715  int num_ents = ents.size();
716  gids = new int[num_ents];
717  ErrorCode result = mbImpl->tag_get_data( mGlobalIdTag, ents, &gids[0] );MB_CHK_SET_ERR( result, "Couldn't get global id data" );
718  minid = *std::min_element( gids, gids + num_ents );
719  maxid = *std::max_element( gids, gids + num_ents );
720  if( 0 == minid )
721  {
722  // gids need to be assigned
723  for( int i = 1; i <= num_ents; i++ )
724  gids[i] = i;
725  result = mbImpl->tag_set_data( mGlobalIdTag, ents, &gids[0] );MB_CHK_SET_ERR( result, "Couldn't set global id data" );
726  maxid = num_ents;
727  }
728 
729  return MB_SUCCESS;
730 }
731 
732 ErrorCode WriteCCMIO::write_nodes( CCMIOID rootID, const Range& verts, const int dimension, CCMIOID& verticesID )
733 {
734  // Get/write map (global ids) first (gids already assigned)
735  unsigned int num_verts = verts.size();
736  std::vector< int > vgids( num_verts );
737  ErrorCode result = mbImpl->tag_get_data( mGlobalIdTag, verts, &vgids[0] );MB_CHK_SET_ERR( result, "Failed to get global ids for vertices" );
738 
739  // Create the map node for vertex ids, and write them to that node
740  CCMIOID mapID;
741  CCMIOError error = kCCMIONoErr;
742  CCMIONewEntity( &error, rootID, kCCMIOMap, "Vertex map", &mapID );
743  CHK_SET_CCMERR( error, "Failure creating Vertex map node" );
744 
745  int maxid = *std::max_element( vgids.begin(), vgids.end() );
746 
747  CCMIOWriteMap( &error, mapID, CCMIOSIZEC( num_verts ), CCMIOSIZEC( maxid ), &vgids[0], CCMIOINDEXC( kCCMIOStart ),
748  CCMIOINDEXC( kCCMIOEnd ) );
749  CHK_SET_CCMERR( error, "Problem writing node map" );
750 
751  // Create the vertex coordinate node, and write it
752  CCMIONewEntity( &error, rootID, kCCMIOVertices, "Vertices", &verticesID );
753  CHK_SET_CCMERR( error, "Trouble creating vertices node" );
754 
755  // Get the vertex locations
756  double* coords = new double[3 * num_verts];
757  std::vector< double* > coord_arrays( 3 );
758  // Cppcheck warning (false positive): variable coord_arrays is assigned a value that is never
759  // used
760  coord_arrays[0] = coords;
761  coord_arrays[1] = coords + num_verts;
762  coord_arrays[2] = ( dimension == 3 ? coords + 2 * num_verts : NULL );
763  result = mWriteIface->get_node_coords( -1, verts.begin(), verts.end(), 3 * num_verts, coords );
764  if( result != MB_SUCCESS )
765  {
766  delete[] coords;
767  return result;
768  }
769 
770  // Transform coordinates, if necessary
771  result = transform_coords( dimension, num_verts, coords );
772  if( result != MB_SUCCESS )
773  {
774  delete[] coords;
775  MB_SET_ERR( result, "Trouble transforming vertex coordinates" );
776  }
777 
778  // Write the vertices
779  CCMIOWriteVerticesd( &error, verticesID, CCMIOSIZEC( dimension ), 1.0, mapID, coords, CCMIOINDEXC( kCCMIOStart ),
780  CCMIOINDEXC( kCCMIOEnd ) );
781  CHK_SET_CCMERR( error, "CCMIOWriteVertices failed" );
782 
783  // Clean up
784  delete[] coords;
785 
786  return MB_SUCCESS;
787 }
788 
789 ErrorCode WriteCCMIO::transform_coords( const int dimension, const int num_nodes, double* coords )
790 {
791  Tag trans_tag;
793  if( result == MB_TAG_NOT_FOUND )
794  return MB_SUCCESS;
795  else if( MB_SUCCESS != result )
796  return result;
797  double trans_matrix[16];
798  const EntityHandle mesh = 0;
799  result = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
800 
801  double* tmp_coords = coords;
802  for( int i = 0; i < num_nodes; i++, tmp_coords += 1 )
803  {
804  double vec1[3] = { 0.0, 0.0, 0.0 };
805  for( int row = 0; row < 3; row++ )
806  {
807  vec1[row] += ( trans_matrix[( row * 4 ) + 0] * coords[0] );
808  vec1[row] += ( trans_matrix[( row * 4 ) + 1] * coords[num_nodes] );
809  if( 3 == dimension ) vec1[row] += ( trans_matrix[( row * 4 ) + 2] * coords[2 * num_nodes] );
810  }
811 
812  coords[0] = vec1[0];
813  coords[num_nodes] = vec1[1];
814  coords[2 * num_nodes] = vec1[2];
815  }
816 
817  return MB_SUCCESS;
818 }
819 
821  std::vector< MaterialSetData >& matset_data,
822  std::vector< NeumannSetData >& neuset_data,
823  Range& /* verts */,
824  CCMIOID& topologyID )
825 {
826  std::vector< int > connect;
827  ErrorCode result;
828  CCMIOID cellMapID, cells;
829  CCMIOError error = kCCMIONoErr;
830 
831  // Don't usually have anywhere near 31 nodes per element
832  connect.reserve( 31 );
834 
835  // Create the topology node, and the cell and cell map nodes
836  CCMIONewEntity( &error, rootID, kCCMIOTopology, "Topology", &topologyID );
837  CHK_SET_CCMERR( error, "Trouble creating topology node" );
838 
839  CCMIONewEntity( &error, rootID, kCCMIOMap, "Cell map", &cellMapID );
840  CHK_SET_CCMERR( error, "Failure creating Cell Map node" );
841 
842  CCMIONewEntity( &error, topologyID, kCCMIOCells, "Cells", &cells );
843  CHK_SET_CCMERR( error, "Trouble creating Cell node under Topology node" );
844 
845  //================================================
846  // Loop over material sets, doing each one at a time
847  //================================================
848  Range all_elems;
849  unsigned int i, num_elems = 0;
850  int max_id = 1;
851  std::vector< int > egids;
852  int tot_elems = 0;
853 
854  for( unsigned int m = 0; m < matset_data.size(); m++ )
855  tot_elems += matset_data[m].elems.size();
856 
857  for( unsigned int m = 0; m < matset_data.size(); m++ )
858  {
859  unsigned int this_num = matset_data[m].elems.size();
860 
861  //================================================
862  // Save all elements being output
863  //================================================
864  all_elems.merge( matset_data[m].elems );
865 
866  //================================================
867  // Assign global ids for elements being written
868  //================================================
869  egids.resize( matset_data[m].elems.size() );
870  for( i = 0; i < this_num; i++ )
871  egids[i] = max_id++;
872  result = mbImpl->tag_set_data( mGlobalIdTag, matset_data[m].elems, &egids[0] );MB_CHK_SET_ERR( result, "Failed to assign global ids for all elements being written" );
873 
874  //================================================
875  // Write cell ids and material types for this matset; reuse egids for cell mat type
876  //================================================
877  CCMIOWriteMap( &error, cellMapID, CCMIOSIZEC( tot_elems ), CCMIOSIZEC( tot_elems ), &egids[0],
878  CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
879  CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
880  CHK_SET_CCMERR( error, "Trouble writing cell map" );
881 
882  if( -1 == matset_data[m].matsetId )
883  {
884  for( i = 0; i < this_num; i++ )
885  egids[i] = m;
886  }
887  else
888  {
889  for( i = 0; i < this_num; i++ )
890  egids[i] = matset_data[m].matsetId;
891  }
892 
893  CCMIOWriteCells( &error, cells, cellMapID, &egids[0], CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
894  CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
895  CHK_SET_CCMERR( error, "Trouble writing Cell node" );
896 
897  //================================================
898  // Write cell entity types
899  //================================================
900  const EntityHandle* conn;
901  int num_conn;
902  int has_mid_nodes[4];
903  std::vector< EntityHandle > storage;
904  for( i = 0, rit = matset_data[m].elems.begin(); i < this_num; i++, ++rit )
905  {
906  result = mbImpl->get_connectivity( *rit, conn, num_conn, false, &storage );MB_CHK_SET_ERR( result, "Trouble getting connectivity for entity type check" );
907  CN::HasMidNodes( mbImpl->type_from_handle( *rit ), num_conn, has_mid_nodes );
908  egids[i] = moab_to_ccmio_type( mbImpl->type_from_handle( *rit ), has_mid_nodes );
909  }
910 
911  CCMIOWriteOpt1i( &error, cells, "CellTopologyType", CCMIOSIZEC( tot_elems ), &egids[0],
912  CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
913  CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
914  CHK_SET_CCMERR( error, "Failed to write cell topo types" );
915 
916  num_elems += this_num;
917  }
918 
919  //================================================
920  // Get skin and neumann set faces
921  //================================================
922  Range neuset_facets, skin_facets;
923  Skinner skinner( mbImpl );
924  result = skinner.find_skin( 0, all_elems, mDimension - 1, skin_facets );MB_CHK_SET_ERR( result, "Failed to get skin facets" );
925 
926  // Remove neumann set facets from skin facets, we have to output these
927  // separately
928  for( i = 0; i < neuset_data.size(); i++ )
929  neuset_facets.merge( neuset_data[i].elems );
930 
931  skin_facets -= neuset_facets;
932  // Make neuset_facets the union, and get ids for them
933  neuset_facets.merge( skin_facets );
934  result = mWriteIface->assign_ids( neuset_facets, mGlobalIdTag, 1 );
935 
936  int fmaxid = neuset_facets.size();
937 
938  //================================================
939  // Write external faces
940  //================================================
941  for( i = 0; i < neuset_data.size(); i++ )
942  {
944  unsigned char cmarks[2];
945  Range ext_faces;
946  std::vector< EntityHandle > mcells;
947  // Removing the faces connected to two regions
948  for( rrit = neuset_data[i].elems.rbegin(); rrit != neuset_data[i].elems.rend(); ++rrit )
949  {
950  mcells.clear();
951  result = mbImpl->get_adjacencies( &( *rrit ), 1, mDimension, false, mcells );MB_CHK_SET_ERR( result, "Trouble getting bounding cells" );
952 
953  result = mbImpl->tag_get_data( mEntityMark, &mcells[0], mcells.size(), cmarks );MB_CHK_SET_ERR( result, "Trouble getting mark tags on cells bounding facets" );
954 
955  if( mcells.size() == 2 && ( mWholeMesh || ( cmarks[0] && cmarks[1] ) ) )
956  {
957  }
958  else
959  {
960  // External face
961  ext_faces.insert( *rrit );
962  }
963  }
964  if( ext_faces.size() != 0 && neuset_data[i].neusetId != 0 )
965  {
966  result = write_external_faces( rootID, topologyID, neuset_data[i].neusetId, ext_faces );MB_CHK_SET_ERR( result, "Trouble writing Neumann set facets" );
967  }
968  ext_faces.clear();
969  }
970 
971  if( !skin_facets.empty() )
972  {
973  result = write_external_faces( rootID, topologyID, 0, skin_facets );MB_CHK_SET_ERR( result, "Trouble writing skin facets" );
974  }
975 
976  //================================================
977  // Now internal faces; loop over elements, do each face on the element
978  //================================================
979  // Mark tag, for face marking on each non-polyhedral element
980 
981  if( num_elems > 1 )
982  { // No internal faces for just one element
983  Tag fmark_tag;
984  unsigned char mval = 0x0, omval;
985  result = mbImpl->tag_get_handle( "__fmark", 1, MB_TYPE_OPAQUE, fmark_tag, MB_TAG_DENSE | MB_TAG_CREAT, &mval );MB_CHK_SET_ERR( result, "Couldn't create mark tag" );
986 
987  std::vector< EntityHandle > tmp_face_cells, storage;
988  std::vector< int > iface_connect, iface_cells;
989  EntityHandle tmp_connect[CN::MAX_NODES_PER_ELEMENT]; // tmp connect vector
990  const EntityHandle *connectc, *oconnectc;
991  int num_connectc; // Cell connectivity
992  const EntityHandle* connectf;
993  int num_connectf; // Face connectivity
994 
995  for( i = 0, rit = all_elems.begin(); i < num_elems; i++, ++rit )
996  {
997  EntityType etype = TYPE_FROM_HANDLE( *rit );
998 
999  //-----------------------
1000  // If not polyh, get mark
1001  //-----------------------
1002  if( MBPOLYHEDRON != etype && MBPOLYGON != etype )
1003  {
1004  result = mbImpl->tag_get_data( fmark_tag, &( *rit ), 1, &mval );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
1005  }
1006 
1007  //-----------------------
1008  // Get cell connectivity, and whether it's a polyhedron
1009  //-----------------------
1010  result = mbImpl->get_connectivity( *rit, connectc, num_connectc, false, &storage );MB_CHK_SET_ERR( result, "Couldn't get entity connectivity" );
1011 
1012  // If polyh, write faces directly
1013  bool is_polyh = ( MBPOLYHEDRON == etype );
1014 
1015  int num_facets = ( is_polyh ? num_connectc : CN::NumSubEntities( etype, mDimension - 1 ) );
1016 
1017  //----------------------------------------------------------
1018  // Loop over each facet of element, outputting it if not marked
1019  //----------------------------------------------------------
1020  for( int f = 0; f < num_facets; f++ )
1021  {
1022  //.............................................
1023  // If this face marked, skip
1024  //.............................................
1025  if( !is_polyh && ( ( mval >> f ) & 0x1 ) ) continue;
1026 
1027  //.................
1028  // Get face connect and adj cells
1029  //.................
1030  if( !is_polyh )
1031  {
1032  // (from CN)
1033  CN::SubEntityConn( connectc, etype, mDimension - 1, f, tmp_connect, num_connectf );
1034  connectf = tmp_connect;
1035  }
1036  else
1037  {
1038  // Directly
1039  result = mbImpl->get_connectivity( connectc[f], connectf, num_connectf, false );MB_CHK_SET_ERR( result, "Couldn't get polyhedron connectivity" );
1040  }
1041 
1042  //............................
1043  // Get adj cells from face connect (same for poly's and not, since both usually
1044  // go through vertices anyway)
1045  //............................
1046  tmp_face_cells.clear();
1047  result = mbImpl->get_adjacencies( connectf, num_connectf, mDimension, false, tmp_face_cells );MB_CHK_SET_ERR( result, "Error getting adj hexes" );
1048 
1049  //...............................
1050  // If this face only bounds one cell, skip, since we exported external faces
1051  // before this loop
1052  //...............................
1053  if( tmp_face_cells.size() != 2 ) continue;
1054 
1055  //.................
1056  // Switch cells so that *rit is always 1st (face connectivity is always written such
1057  // that that one is with forward sense)
1058  //.................
1059  int side_num = 0, sense = 0, offset = 0;
1060  if( !is_polyh && tmp_face_cells[0] != *rit )
1061  {
1062  EntityHandle tmph = tmp_face_cells[0];
1063  tmp_face_cells[0] = tmp_face_cells[1];
1064  tmp_face_cells[1] = tmph;
1065  }
1066 
1067  //.................
1068  // Save ids of cells
1069  //.................
1070  assert( tmp_face_cells[0] != tmp_face_cells[1] );
1071  iface_cells.resize( iface_cells.size() + 2 );
1072  result = mbImpl->tag_get_data( mGlobalIdTag, &tmp_face_cells[0], tmp_face_cells.size(),
1073  &iface_cells[iface_cells.size() - 2] );MB_CHK_SET_ERR( result, "Trouble getting global ids for bounded cells" );
1074  iface_connect.push_back( num_connectf );
1075 
1076  //.................
1077  // Save indices of face vertices
1078  //.................
1079  unsigned int tmp_size = iface_connect.size();
1080  iface_connect.resize( tmp_size + num_connectf );
1081  result = mbImpl->tag_get_data( mGlobalIdTag, connectf, num_connectf, &iface_connect[tmp_size] );MB_CHK_SET_ERR( result, "Trouble getting global id for internal face" );
1082 
1083  //.................
1084  // Mark other cell with the right side #
1085  //.................
1086  if( !is_polyh )
1087  {
1088  // Mark other cell for this face, if there is another cell
1089 
1090  result = mbImpl->get_connectivity( tmp_face_cells[1], oconnectc, num_connectc, false, &storage );MB_CHK_SET_ERR( result, "Couldn't get other entity connectivity" );
1091 
1092  // Get side number in other cell
1093  CN::SideNumber( TYPE_FROM_HANDLE( tmp_face_cells[1] ), oconnectc, connectf, num_connectf,
1094  mDimension - 1, side_num, sense, offset );
1095  // Set mark for that face on the other cell
1096  result = mbImpl->tag_get_data( fmark_tag, &tmp_face_cells[1], 1, &omval );MB_CHK_SET_ERR( result, "Couldn't get mark data for other cell" );
1097  }
1098 
1099  omval |= ( 0x1 << (unsigned int)side_num );
1100  result = mbImpl->tag_set_data( fmark_tag, &tmp_face_cells[1], 1, &omval );MB_CHK_SET_ERR( result, "Couldn't set mark data for other cell" );
1101  } // Loop over faces in elem
1102  } // Loop over elems
1103 
1104  //================================================
1105  // Write internal faces
1106  //================================================
1107 
1108  CCMIOID mapID;
1109  CCMIONewEntity( &error, rootID, kCCMIOMap, NULL, &mapID );
1110  CHK_SET_CCMERR( error, "Trouble creating Internal Face map node" );
1111 
1112  unsigned int num_ifaces = iface_cells.size() / 2;
1113 
1114  // Set gids for internal faces; reuse egids
1115  egids.resize( num_ifaces );
1116  for( i = 1; i <= num_ifaces; i++ )
1117  egids[i - 1] = fmaxid + i;
1118  CCMIOWriteMap( &error, mapID, CCMIOSIZEC( num_ifaces ), CCMIOSIZEC( fmaxid + num_ifaces ), &egids[0],
1119  CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
1120  CHK_SET_CCMERR( error, "Trouble writing Internal Face map node" );
1121 
1122  CCMIOID id;
1123  CCMIONewEntity( &error, topologyID, kCCMIOInternalFaces, "Internal faces", &id );
1124  CHK_SET_CCMERR( error, "Failed to create Internal face node under Topology node" );
1125  CCMIOWriteFaces( &error, id, kCCMIOInternalFaces, mapID, CCMIOSIZEC( iface_connect.size() ), &iface_connect[0],
1126  CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
1127  CHK_SET_CCMERR( error, "Failure writing Internal face connectivity" );
1128  CCMIOWriteFaceCells( &error, id, kCCMIOInternalFaces, mapID, &iface_cells[0], CCMIOINDEXC( kCCMIOStart ),
1129  CCMIOINDEXC( kCCMIOEnd ) );
1130  CHK_SET_CCMERR( error, "Failure writing Internal face cells" );
1131  }
1132 
1133  return MB_SUCCESS;
1134 }
1135 
1136 int WriteCCMIO::moab_to_ccmio_type( EntityType etype, int has_mid_nodes[] )
1137 {
1138  int ctype = -1;
1139  if( has_mid_nodes[0] || has_mid_nodes[2] || has_mid_nodes[3] ) return ctype;
1140 
1141  switch( etype )
1142  {
1143  case MBVERTEX:
1144  ctype = 1;
1145  break;
1146  case MBEDGE:
1147  if( !has_mid_nodes[1] )
1148  ctype = 2;
1149  else
1150  ctype = 28;
1151  break;
1152  case MBQUAD:
1153  if( has_mid_nodes[1] )
1154  ctype = 4;
1155  else
1156  ctype = 3;
1157  break;
1158  case MBTET:
1159  if( has_mid_nodes[1] )
1160  ctype = 23;
1161  else
1162  ctype = 13;
1163  break;
1164  case MBPRISM:
1165  if( has_mid_nodes[1] )
1166  ctype = 22;
1167  else
1168  ctype = 12;
1169  break;
1170  case MBPYRAMID:
1171  if( has_mid_nodes[1] )
1172  ctype = 24;
1173  else
1174  ctype = 14;
1175  break;
1176  case MBHEX:
1177  if( has_mid_nodes[1] )
1178  ctype = 21;
1179  else
1180  ctype = 11;
1181  break;
1182  case MBPOLYHEDRON:
1183  ctype = 255;
1184  break;
1185  default:
1186  break;
1187  }
1188 
1189  return ctype;
1190 }
1191 
1192 ErrorCode WriteCCMIO::write_external_faces( CCMIOID rootID, CCMIOID topologyID, int set_num, Range& facets )
1193 {
1194  CCMIOError error = kCCMIONoErr;
1195  CCMIOID mapID, id;
1196 
1197  // Get gids for these faces
1198  int *gids = NULL, minid, maxid;
1199  ErrorCode result = get_gids( facets, gids, minid, maxid );MB_CHK_SET_ERR( result, "Trouble getting global ids for facets" );
1200 
1201  // Write the face id map
1202  CCMIONewEntity( &error, rootID, kCCMIOMap, NULL, &mapID );
1203  CHK_SET_CCMERR( error, "Problem creating face id map" );
1204 
1205  CCMIOWriteMap( &error, mapID, CCMIOSIZEC( facets.size() ), CCMIOSIZEC( maxid ), gids, CCMIOINDEXC( kCCMIOStart ),
1206  CCMIOINDEXC( kCCMIOEnd ) );
1207  CHK_SET_CCMERR( error, "Problem writing face id map" );
1208 
1209  // Get the connectivity of the faces; set size by how many verts in last facet
1210  const EntityHandle* connect;
1211  int num_connect;
1212  result = mbImpl->get_connectivity( *facets.rbegin(), connect, num_connect );MB_CHK_SET_ERR( result, "Failed to get connectivity of last facet" );
1213  std::vector< int > fconnect( facets.size() * ( num_connect + 1 ) );
1214 
1215  result = mWriteIface->get_element_connect( facets.begin(), facets.end(), num_connect, mGlobalIdTag, fconnect.size(),
1216  &fconnect[0], true );MB_CHK_SET_ERR( result, "Failed to get facet connectivity" );
1217 
1218  // Get and write a new external face entity
1219  CCMIONewIndexedEntity( &error, topologyID, kCCMIOBoundaryFaces, set_num, "Boundary faces", &id );
1220  CHK_SET_CCMERR( error, "Problem creating boundary face entity" );
1221 
1222  CCMIOWriteFaces( &error, id, kCCMIOBoundaryFaces, mapID, CCMIOSIZEC( fconnect.size() ), &fconnect[0],
1223  CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
1224  CHK_SET_CCMERR( error, "Problem writing boundary faces" );
1225 
1226  // Get info on bounding cells; reuse fconnect
1227  std::vector< EntityHandle > cells;
1228  unsigned char cmarks[2];
1229  int i, j = 0;
1230  Range dead_facets;
1231  Range::iterator rit;
1232 
1233  // About error checking in this loop: if any facets have no bounding cells,
1234  // this is an error, since global ids for facets are computed outside this loop
1235  for( rit = facets.begin(), i = 0; rit != facets.end(); ++rit, i++ )
1236  {
1237  cells.clear();
1238 
1239  // Get cell then gid of cell
1240  result = mbImpl->get_adjacencies( &( *rit ), 1, mDimension, false, cells );MB_CHK_SET_ERR( result, "Trouble getting bounding cells" );
1241  if( cells.empty() )
1242  {
1243  MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with no output bounding cell" );
1244  }
1245 
1246  // Check we don't bound more than one cell being output
1247  result = mbImpl->tag_get_data( mEntityMark, &cells[0], cells.size(), cmarks );MB_CHK_SET_ERR( result, "Trouble getting mark tags on cells bounding facets" );
1248  if( cells.size() == 2 && ( mWholeMesh || ( cmarks[0] && cmarks[1] ) ) )
1249  {
1250  MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with two output bounding cells" );
1251  }
1252  else if( 1 == cells.size() && !mWholeMesh && !cmarks[0] )
1253  {
1254  MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with no output bounding cells" );
1255  }
1256 
1257  // Make sure 1st cell is the one being output
1258  if( 2 == cells.size() && !( cmarks[0] | 0x0 ) && ( cmarks[1] & 0x1 ) ) cells[0] = cells[1];
1259 
1260  // Get gid for bounded cell
1261  result = mbImpl->tag_get_data( mGlobalIdTag, &cells[0], 1, &fconnect[j] );MB_CHK_SET_ERR( result, "Couldn't get global id tag for bounded cell" );
1262 
1263  j++;
1264  }
1265 
1266  // Write the bounding cell data
1267  CCMIOWriteFaceCells( &error, id, kCCMIOBoundaryFaces, mapID, &fconnect[0], CCMIOINDEXC( kCCMIOStart ),
1268  CCMIOINDEXC( kCCMIOEnd ) );
1269  CHK_SET_CCMERR( error, "Problem writing boundary cell data" );
1270 
1271  return MB_SUCCESS;
1272 }
1273 
1275  int current_sense,
1276  Range& forward_elems,
1277  Range& reverse_elems )
1278 {
1279  Range neuset_elems, neuset_meshsets;
1280 
1281  // Get the sense tag; don't need to check return, might be an error if the tag
1282  // hasn't been created yet
1283  Tag sense_tag = 0;
1284  mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
1285 
1286  // Get the entities in this set, non-recursive
1287  ErrorCode result = mbImpl->get_entities_by_handle( neuset, neuset_elems );
1288  if( MB_FAILURE == result ) return result;
1289 
1290  // Now remove the meshsets into the neuset_meshsets; first find the first meshset,
1291  Range::iterator range_iter = neuset_elems.begin();
1292  while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != neuset_elems.end() )
1293  ++range_iter;
1294 
1295  // Then, if there are some, copy them into neuset_meshsets and erase from neuset_elems
1296  if( range_iter != neuset_elems.end() )
1297  {
1298  std::copy( range_iter, neuset_elems.end(), range_inserter( neuset_meshsets ) );
1299  neuset_elems.erase( range_iter, neuset_elems.end() );
1300  }
1301 
1302  // OK, for the elements, check the sense of this set and copy into the right range
1303  // (if the sense is 0, copy into both ranges)
1304 
1305  // Need to step forward on list until we reach the right dimension
1306  Range::iterator dum_it = neuset_elems.end();
1307  --dum_it;
1308  int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
1309  dum_it = neuset_elems.begin();
1310  while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != neuset_elems.end() )
1311  ++dum_it;
1312 
1313  if( current_sense == 1 || current_sense == 0 )
1314  std::copy( dum_it, neuset_elems.end(), range_inserter( forward_elems ) );
1315  if( current_sense == -1 || current_sense == 0 )
1316  std::copy( dum_it, neuset_elems.end(), range_inserter( reverse_elems ) );
1317 
1318  // Now loop over the contained meshsets, getting the sense of those and calling this
1319  // function recursively
1320  for( range_iter = neuset_meshsets.begin(); range_iter != neuset_meshsets.end(); ++range_iter )
1321  {
1322  // First get the sense; if it's not there, by convention it's forward
1323  int this_sense;
1324  if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
1325  this_sense = 1;
1326 
1327  // Now get all the entities on this meshset, with the proper (possibly reversed) sense
1328  get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
1329  }
1330 
1331  return result;
1332 }
1333 } // namespace moab