Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadCCMIO.cpp
Go to the documentation of this file.
1 #include <cstdlib> // For exit()
2 #include <vector>
3 #include <map>
4 #include <iostream>
5 #include <string>
6 #include <algorithm>
7 
8 #include "moab/CN.hpp"
9 #include "moab/Range.hpp"
10 #include "moab/Interface.hpp"
11 #include "MBTagConventions.hpp"
12 #include "Internals.hpp"
13 #include "moab/ReadUtilIface.hpp"
14 #include "moab/FileOptions.hpp"
15 #include "ReadCCMIO.hpp"
16 #include "moab/MeshTopoUtil.hpp"
17 
18 #include "ccmio.h"
19 
20 /*
21  * CCMIO file structure
22  *
23  * Root
24  * State(kCCMIOState)
25  * Processor*
26  * Vertices
27  * ->ReadVerticesx, ReadMap
28  * Topology
29  * Boundary faces*(kCCMIOBoundaryFaces)
30  * ->ReadFaces, ReadFaceCells, ReadMap
31  * Internal faces(kCCMIOInternalFaces)
32  * Cells (kCCMIOCells)
33  * ->ReadCells (mapID), ReadMap, ReadCells (cellTypes)
34  * Solution
35  * Phase
36  * Field
37  * FieldData
38  * Problem(kCCMIOProblemDescription)
39  * CellType* (kCCMIOCellType)
40  * Index (GetEntityIndex), MaterialId(ReadOpti), MaterialType(ReadOptstr),
41  * PorosityId(ReadOpti), SpinId(ReadOpti), GroupId(ReadOpti)
42  *
43  * MaterialType (CCMIOReadOptstr in readexample)
44  * constants (see readexample)
45  * lagrangian data (CCMIOReadLagrangianData)
46  * vertices label (CCMIOEntityDescription)
47  * restart info: char solver[], iteratoins, time, char timeUnits[], angle
48  * (CCMIOReadRestartInfo, kCCMIORestartData), reference data?
49  * phase:
50  * field: char name[], dims, CCMIODataType datatype, char units[]
51  * dims = kCCMIOScalar (CCMIOReadFieldDataf),
52  * kCCMIOVector (CCMIOReadMultiDimensionalFieldData),
53  * kCCMIOTensor
54  * MonitoringSets: num, name (CellSet, VertexSet, BoundarySet, BlockSet, SplineSet, CoupleSet)
55  * CCMIOGetProstarSet, CCMIOReadOpt1i,
56  */
57 
58 enum DataType
59 {
60  kScalar,
61  kVector,
62  kVertex,
63  kCell,
64  kInternalFace,
65  kBoundaryFace,
66  kBoundaryData,
67  kBoundaryFaceData,
68  kCellType
69 };
70 
71 namespace moab
72 {
73 
74 static int const kNValues = 10; // Number of values of each element to print
75 static char const kDefaultState[] = "default";
76 static char const kUnitsName[] = "Units";
77 static int const kVertOffset = 2;
78 static int const kCellInc = 4;
79 
80 #define CHK_SET_CCMERR( ccm_err_code, ccm_err_msg ) \
81  { \
82  if( kCCMIONoErr != ( ccm_err_code ) && kCCMIONoFileErr != ( ccm_err_code ) && \
83  kCCMIONoNodeErr != ( ccm_err_code ) ) \
84  MB_SET_ERR( MB_FAILURE, ccm_err_msg ); \
85  }
86 
87 ReaderIface* ReadCCMIO::factory( Interface* iface )
88 {
89  return new ReadCCMIO( iface );
90 }
91 
92 ReadCCMIO::ReadCCMIO( Interface* impl )
93  : mMaterialIdTag( 0 ), mMaterialTypeTag( 0 ), mRadiationTag( 0 ), mPorosityIdTag( 0 ), mSpinIdTag( 0 ),
94  mGroupIdTag( 0 ), mColorIdxTag( 0 ), mProcessorIdTag( 0 ), mLightMaterialTag( 0 ), mFreeSurfaceMaterialTag( 0 ),
95  mThicknessTag( 0 ), mProstarRegionNumberTag( 0 ), mBoundaryTypeTag( 0 ), mCreatingProgramTag( 0 ), mbImpl( impl ),
96  hasSolution( false )
97 {
98  assert( impl != NULL );
99 
100  impl->query_interface( readMeshIface );
101 
102  // Initialize in case tag_get_handle fails below
103  mMaterialSetTag = 0;
104  mDirichletSetTag = 0;
105  mNeumannSetTag = 0;
106  mHasMidNodesTag = 0;
107  mGlobalIdTag = 0;
108  mNameTag = 0;
109 
110  //! Get and cache predefined tag handles
111  const int negone = -1;
112  ErrorCode result = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
113  MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get MATERIAL_SET tag" );
114 
115  result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
116  MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get DIRICHLET_SET tag" );
117 
118  result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
119  MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get NEUMANN_SET tag" );
120 
121  const int negonearr[] = { -1, -1, -1, -1 };
122  result = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
123  MB_TAG_CREAT | MB_TAG_SPARSE, negonearr );MB_CHK_SET_ERR_RET( result, "Failed to get HAS_MID_NODES tag" );
124 
125  mGlobalIdTag = impl->globalId_tag();
126 
127  result =
128  impl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, mNameTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_RET( result, "Failed to get NAME tag" );
129 }
130 
131 ReadCCMIO::~ReadCCMIO()
132 {
133  mbImpl->release_interface( readMeshIface );
134 }
135 
136 ErrorCode ReadCCMIO::load_file( const char* file_name,
137  const EntityHandle* file_set,
138  const FileOptions& /* opts */,
139  const ReaderIface::SubsetList* subset_list,
140  const Tag* /* file_id_tag */ )
141 {
142  CCMIOID rootID, problemID, stateID, processorID, verticesID, topologyID, solutionID;
143  CCMIOError error = kCCMIONoErr;
144 
145  if( subset_list )
146  {
147  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for CCMOI data" );
148  }
149 
150  CCMIOOpenFile( &error, file_name, kCCMIORead, &rootID );
151  CHK_SET_CCMERR( error, "Problem opening file" );
152 
153  // Get the file state
154  ErrorCode rval = get_state( rootID, problemID, stateID );MB_CHK_SET_ERR( rval, "Failed to get state" );
155 
156  // Get processors
157  std::vector< CCMIOSize_t > procs;
158  bool has_solution = false;
159  rval = get_processors( stateID, processorID, verticesID, topologyID, solutionID, procs, has_solution );MB_CHK_SET_ERR( rval, "Failed to get processors" );
160 
161  std::vector< CCMIOSize_t >::iterator vit;
162  Range new_ents, *new_ents_ptr = NULL;
163  if( file_set ) new_ents_ptr = &new_ents;
164 
165  for( vit = procs.begin(); vit != procs.end(); ++vit )
166  {
167  rval = read_processor( stateID, problemID, processorID, verticesID, topologyID, *vit, new_ents_ptr );MB_CHK_SET_ERR( rval, "Failed to read processors" );
168  }
169 
170  // Load some meta-data
171  rval = load_metadata( rootID, problemID, stateID, processorID, file_set );MB_CHK_SET_ERR( rval, "Failed to load some meta-data" );
172 
173  // Now, put all this into the file set, if there is one
174  if( file_set )
175  {
176  rval = mbImpl->add_entities( *file_set, new_ents );MB_CHK_SET_ERR( rval, "Failed to add new entities to file set" );
177  }
178 
179  return rval;
180 }
181 
182 ErrorCode ReadCCMIO::get_state( CCMIOID rootID, CCMIOID& problemID, CCMIOID& stateID )
183 {
184  CCMIOError error = kCCMIONoErr;
185 
186  // First try default
187  CCMIOGetState( &error, rootID, "default", &problemID, &stateID );
188  if( kCCMIONoErr != error )
189  {
190  CCMIOSize_t i = CCMIOSIZEC( 0 );
191  CCMIOError tmp_error = kCCMIONoErr;
192  CCMIONextEntity( &tmp_error, rootID, kCCMIOState, &i, &stateID );
193  if( kCCMIONoErr == tmp_error ) CCMIONextEntity( &error, rootID, kCCMIOProblemDescription, &i, &problemID );
194  }
195  CHK_SET_CCMERR( error, "Couldn't find state" );
196 
197  return MB_SUCCESS;
198 }
199 
200 ErrorCode ReadCCMIO::load_metadata( CCMIOID rootID,
201  CCMIOID problemID,
202  CCMIOID /* stateID */,
203  CCMIOID processorID,
204  const EntityHandle* file_set )
205 {
206  // Read the simulation title.
207  CCMIOError error = kCCMIONoErr;
208  ErrorCode rval = MB_SUCCESS;
209  CCMIONode rootNode;
210  if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
211  {
212  char* name = NULL;
213  CCMIOGetTitle( &error, rootNode, &name );
214 
215  if( NULL != name && strlen( name ) != 0 )
216  {
217  // Make a tag for it and tag the read set
218  Tag simname;
219  rval = mbImpl->tag_get_handle( "Title", strlen( name ), MB_TYPE_OPAQUE, simname,
220  MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Simulation name tag not found or created" );
221  EntityHandle set = file_set ? *file_set : 0;
222  rval = mbImpl->tag_set_data( simname, &set, 1, name );MB_CHK_SET_ERR( rval, "Problem setting simulation name tag" );
223  }
224  if( name ) free( name );
225  }
226 
227  // Creating program
228  EntityHandle dumh = ( file_set ? *file_set : 0 );
229  rval = get_str_option( "CreatingProgram", dumh, mCreatingProgramTag, processorID );MB_CHK_SET_ERR( rval, "Trouble getting CreatingProgram tag" );
230 
231  rval = load_matset_data( problemID );MB_CHK_SET_ERR( rval, "Failure loading matset data" );
232 
233  rval = load_neuset_data( problemID );MB_CHK_SET_ERR( rval, "Failure loading neuset data" );
234 
235  return rval;
236 }
237 
238 ErrorCode ReadCCMIO::load_matset_data( CCMIOID problemID )
239 {
240  // Make sure there are matsets
241  if( newMatsets.empty() ) return MB_SUCCESS;
242 
243  // ... walk through each cell type
244  CCMIOSize_t i = CCMIOSIZEC( 0 );
245  CCMIOID next;
246  CCMIOError error = kCCMIONoErr;
247 
248  while( CCMIONextEntity( NULL, problemID, kCCMIOCellType, &i, &next ) == kCCMIONoErr )
249  {
250  // Get index, corresponding set, and label with material set tag
251  int mindex;
252  CCMIOGetEntityIndex( &error, next, &mindex );
253  std::map< int, EntityHandle >::iterator mit = newMatsets.find( mindex );
254  if( mit == newMatsets.end() )
255  // No actual faces for this matset; continue to next
256  continue;
257 
258  EntityHandle dum_ent = mit->second;
259  ErrorCode rval = mbImpl->tag_set_data( mMaterialSetTag, &dum_ent, 1, &mindex );MB_CHK_SET_ERR( rval, "Trouble setting material set tag" );
260 
261  // Set name
262  CCMIOSize_t len;
263  CCMIOEntityLabel( &error, next, &len, NULL );
264  std::vector< char > opt_string2( GETINT32( len ) + 1, '\0' );
265  CCMIOEntityLabel( &error, next, NULL, &opt_string2[0] );
266  if( opt_string2.size() >= NAME_TAG_SIZE )
267  opt_string2[NAME_TAG_SIZE - 1] = '\0';
268  else
269  ( opt_string2.resize( NAME_TAG_SIZE, '\0' ) );
270  rval = mbImpl->tag_set_data( mNameTag, &dum_ent, 1, &opt_string2[0] );MB_CHK_SET_ERR( rval, "Trouble setting name tag for material set" );
271 
272  // Material id
273  rval = get_int_option( "MaterialId", dum_ent, mMaterialIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting MaterialId tag" );
274 
275  rval = get_str_option( "MaterialType", dum_ent, mMaterialTypeTag, next );MB_CHK_SET_ERR( rval, "Trouble getting MaterialType tag" );
276 
277  rval = get_int_option( "Radiation", dum_ent, mRadiationTag, next );MB_CHK_SET_ERR( rval, "Trouble getting Radiation option" );
278 
279  rval = get_int_option( "PorosityId", dum_ent, mPorosityIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting PorosityId option" );
280 
281  rval = get_int_option( "SpinId", dum_ent, mSpinIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting SpinId option" );
282 
283  rval = get_int_option( "GroupId", dum_ent, mGroupIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting GroupId option" );
284 
285  rval = get_int_option( "ColorIdx", dum_ent, mColorIdxTag, next );MB_CHK_SET_ERR( rval, "Trouble getting ColorIdx option" );
286 
287  rval = get_int_option( "ProcessorId", dum_ent, mProcessorIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting ProcessorId option" );
288 
289  rval = get_int_option( "LightMaterial", dum_ent, mLightMaterialTag, next );MB_CHK_SET_ERR( rval, "Trouble getting LightMaterial option" );
290 
291  rval = get_int_option( "FreeSurfaceMaterial", dum_ent, mFreeSurfaceMaterialTag, next );MB_CHK_SET_ERR( rval, "Trouble getting FreeSurfaceMaterial option" );
292 
293  rval = get_dbl_option( "Thickness", dum_ent, mThicknessTag, next );MB_CHK_SET_ERR( rval, "Trouble getting Thickness option" );
294  }
295 
296  return MB_SUCCESS;
297 }
298 
299 ErrorCode ReadCCMIO::get_int_option( const char* opt_str, EntityHandle seth, Tag& tag, CCMIOID node )
300 {
301  int idum;
302  ErrorCode rval;
303  if( kCCMIONoErr == CCMIOReadOpti( NULL, node, opt_str, &idum ) )
304  {
305  if( !tag )
306  {
307  rval = mbImpl->tag_get_handle( opt_str, 1, MB_TYPE_INTEGER, tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
308  }
309 
310  rval = mbImpl->tag_set_data( tag, &seth, 1, &idum );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
311  }
312 
313  return MB_SUCCESS;
314 }
315 
316 ErrorCode ReadCCMIO::get_dbl_option( const char* opt_str, EntityHandle seth, Tag& tag, CCMIOID node )
317 {
318  float fdum;
319  if( kCCMIONoErr == CCMIOReadOptf( NULL, node, opt_str, &fdum ) )
320  {
321  ErrorCode rval;
322  if( !tag )
323  {
324  rval = mbImpl->tag_get_handle( opt_str, 1, MB_TYPE_DOUBLE, tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
325  }
326 
327  double dum_dbl = fdum;
328  rval = mbImpl->tag_set_data( tag, &seth, 1, &dum_dbl );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
329  }
330 
331  return MB_SUCCESS;
332 }
333 
334 ErrorCode ReadCCMIO::get_str_option( const char* opt_str,
335  EntityHandle seth,
336  Tag& tag,
337  CCMIOID node,
338  const char* other_tag_name )
339 {
340  int len;
341  CCMIOError error = kCCMIONoErr;
342  std::vector< char > opt_string;
343  if( kCCMIONoErr != CCMIOReadOptstr( NULL, node, opt_str, &len, NULL ) ) return MB_SUCCESS;
344 
345  opt_string.resize( len );
346  CCMIOReadOptstr( &error, node, opt_str, &len, &opt_string[0] );
347  ErrorCode rval = MB_SUCCESS;
348  if( !tag )
349  {
350  rval = mbImpl->tag_get_handle( other_tag_name ? other_tag_name : opt_str, NAME_TAG_SIZE, MB_TYPE_OPAQUE, tag,
351  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
352  }
353 
354  if( opt_string.size() > NAME_TAG_SIZE )
355  opt_string[NAME_TAG_SIZE - 1] = '\0';
356  else
357  ( opt_string.resize( NAME_TAG_SIZE, '\0' ) );
358 
359  rval = mbImpl->tag_set_data( tag, &seth, 1, &opt_string[0] );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
360 
361  return MB_SUCCESS;
362 }
363 
364 ErrorCode ReadCCMIO::load_neuset_data( CCMIOID problemID )
365 {
366  CCMIOSize_t i = CCMIOSIZEC( 0 );
367  CCMIOID next;
368 
369  // Make sure there are matsets
370  if( newNeusets.empty() ) return MB_SUCCESS;
371 
372  while( CCMIONextEntity( NULL, problemID, kCCMIOBoundaryRegion, &i, &next ) == kCCMIONoErr )
373  {
374  // Get index, corresponding set, and label with neumann set tag
375  int mindex;
376  CCMIOError error = kCCMIONoErr;
377  CCMIOGetEntityIndex( &error, next, &mindex );
378  std::map< int, EntityHandle >::iterator mit = newNeusets.find( mindex );
379  if( mit == newNeusets.end() )
380  // No actual faces for this neuset; continue to next
381  continue;
382 
383  EntityHandle dum_ent = mit->second;
384  ErrorCode rval = mbImpl->tag_set_data( mNeumannSetTag, &dum_ent, 1, &mindex );MB_CHK_SET_ERR( rval, "Trouble setting neumann set tag" );
385 
386  // Set name
387  rval = get_str_option( "BoundaryName", dum_ent, mNameTag, next, NAME_TAG_NAME );MB_CHK_SET_ERR( rval, "Trouble creating BoundaryName tag" );
388 
389  // BoundaryType
390  rval = get_str_option( "BoundaryType", dum_ent, mBoundaryTypeTag, next );MB_CHK_SET_ERR( rval, "Trouble creating BoundaryType tag" );
391 
392  // ProstarRegionNumber
393  rval = get_int_option( "ProstarRegionNumber", dum_ent, mProstarRegionNumberTag, next );MB_CHK_SET_ERR( rval, "Trouble creating ProstarRegionNumber tag" );
394  }
395 
396  return MB_SUCCESS;
397 }
398 
399 ErrorCode ReadCCMIO::read_processor( CCMIOID /* stateID */,
400  CCMIOID problemID,
401  CCMIOID processorID,
402  CCMIOID verticesID,
403  CCMIOID topologyID,
404  CCMIOSize_t proc,
405  Range* new_ents )
406 {
407  ErrorCode rval;
408 
409  // vert_map fields: s: none, i: gid, ul: vert handle, r: none
410  // TupleList vert_map(0, 1, 1, 0, 0);
411  TupleList vert_map;
412  rval = read_vertices( proc, processorID, verticesID, topologyID, new_ents, vert_map );MB_CHK_SET_ERR( rval, "Failed to read vertices" );
413 
414  rval = read_cells( proc, problemID, verticesID, topologyID, vert_map, new_ents );MB_CHK_SET_ERR( rval, "Failed to read cells" );
415 
416  return rval;
417 }
418 
419 ErrorCode ReadCCMIO::read_cells( CCMIOSize_t /* proc */,
420  CCMIOID problemID,
421  CCMIOID /* verticesID */,
422  CCMIOID topologyID,
423  TupleList& vert_map,
424  Range* new_ents )
425 {
426  // Read the faces.
427  // face_map fields: s:forward/reverse, i: cell id, ul: face handle, r: none
428  ErrorCode rval;
429 #ifdef TUPLE_LIST
430  TupleList face_map( 1, 1, 1, 0, 0 );
431 #else
432  TupleList face_map;
433  SenseList sense_map;
434 #endif
435  rval = read_all_faces( topologyID, vert_map, face_map,
436 #ifndef TUPLE_LIST
437  sense_map,
438 #endif
439  new_ents );MB_CHK_SET_ERR( rval, "Failed to read all cells" );
440 
441  // Read the cell topology types, if any exist in the file
442  std::map< int, int > cell_topo_types;
443  rval = read_topology_types( topologyID, cell_topo_types );MB_CHK_SET_ERR( rval, "Problem reading cell topo types" );
444 
445  // Now construct the cells; sort the face map by cell ids first
446 #ifdef TUPLE_LIST
447  rval = face_map.sort( 1 );MB_CHK_SET_ERR( rval, "Couldn't sort face map by cell id" );
448 #endif
449  std::vector< EntityHandle > new_cells;
450  rval = construct_cells( face_map,
451 #ifndef TUPLE_LIST
452  sense_map,
453 #endif
454  vert_map, cell_topo_types, new_cells );MB_CHK_SET_ERR( rval, "Failed to construct cells" );
455  if( new_ents )
456  {
457  Range::iterator rit = new_ents->end();
458  std::vector< EntityHandle >::reverse_iterator vit;
459  for( vit = new_cells.rbegin(); vit != new_cells.rend(); ++vit )
460  rit = new_ents->insert( rit, *vit );
461  }
462 
463  rval = read_gids_and_types( problemID, topologyID, new_cells );MB_CHK_SET_ERR( rval, "Failed to read gids and types" );
464 
465  return MB_SUCCESS;
466 }
467 
468 ErrorCode ReadCCMIO::read_topology_types( CCMIOID& topologyID, std::map< int, int >& cell_topo_types )
469 {
470  CCMIOError error = kCCMIONoErr;
471  CCMIOID cellID, mapID;
472  CCMIOSize_t ncells;
473  CCMIOGetEntity( &error, topologyID, kCCMIOCells, 0, &cellID );
474  CCMIOEntitySize( &error, cellID, &ncells, NULL );
475  int num_cells = GETINT32( ncells );
476 
477  // First, do a dummy read to see if we even have topo types in this mesh
478  int dum_int;
479  CCMIOReadOpt1i( &error, cellID, "CellTopologyType", &dum_int, CCMIOINDEXC( kCCMIOStart ),
480  CCMIOINDEXC( kCCMIOStart ) + 1 );
481  if( kCCMIONoErr != error ) return MB_SUCCESS;
482 
483  // OK, we have topo types; first get the map node
484  std::vector< int > dum_ints( num_cells );
485  CCMIOReadCells( &error, cellID, &mapID, &dum_ints[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOStart ) + 1 );
486  CHK_SET_CCMERR( error, "Failed to get the map node" );
487 
488  // Now read the map
489  CCMIOReadMap( &error, mapID, &dum_ints[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
490  CHK_SET_CCMERR( error, "Failed to get cell ids" );
491  int i;
492  for( i = 0; i < num_cells; i++ )
493  cell_topo_types[dum_ints[i]] = 0;
494 
495  // Now read the cell topo types for real, reusing cell_topo_types
496  std::vector< int > topo_types( num_cells );
497  CCMIOReadOpt1i( &error, cellID, "CellTopologyType", &topo_types[0], CCMIOINDEXC( kCCMIOStart ),
498  CCMIOINDEXC( kCCMIOEnd ) );
499  CHK_SET_CCMERR( error, "Failed to get cell topo types" );
500  for( i = 0; i < num_cells; i++ )
501  cell_topo_types[dum_ints[i]] = topo_types[i];
502 
503  return MB_SUCCESS;
504 }
505 
506 ErrorCode ReadCCMIO::read_gids_and_types( CCMIOID /* problemID */,
507  CCMIOID topologyID,
508  std::vector< EntityHandle >& cells )
509 {
510  // Get the cells entity and number of cells
511  CCMIOSize_t dum_cells;
512  int num_cells;
513  CCMIOError error = kCCMIONoErr;
514  CCMIOID cellsID, mapID;
515  CCMIOGetEntity( &error, topologyID, kCCMIOCells, 0, &cellsID );
516  CCMIOEntitySize( &error, cellsID, &dum_cells, NULL );
517  num_cells = GETINT32( dum_cells );
518 
519  // Check the number of cells against how many are in the cell array
520  if( num_cells != (int)cells.size() ) MB_SET_ERR( MB_FAILURE, "Number of cells doesn't agree" );
521 
522  // Read the gid map and set global ids
523  std::vector< int > cell_gids( num_cells );
524  CCMIOReadCells( &error, cellsID, &mapID, NULL, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
525  CHK_SET_CCMERR( error, "Couldn't read cells" );
526  CCMIOReadMap( &error, mapID, &cell_gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
527  CHK_SET_CCMERR( error, "Couldn't read cell id map" );
528 
529  ErrorCode rval = mbImpl->tag_set_data( mGlobalIdTag, &cells[0], cells.size(), &cell_gids[0] );MB_CHK_SET_ERR( rval, "Couldn't set gids tag" );
530 
531  // Now read cell material types; reuse cell_gids
532  CCMIOReadCells( &error, cellsID, NULL, &cell_gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
533  CHK_SET_CCMERR( error, "Trouble reading cell types" );
534 
535  // Create the matsets
536  std::map< int, Range > matset_ents;
537  for( int i = 0; i < num_cells; i++ )
538  matset_ents[cell_gids[i]].insert( cells[i] );
539 
540  for( std::map< int, Range >::iterator mit = matset_ents.begin(); mit != matset_ents.end(); ++mit )
541  {
542  EntityHandle matset;
543  rval = mbImpl->create_meshset( MESHSET_SET, matset );MB_CHK_SET_ERR( rval, "Couldn't create material set" );
544  newMatsets[mit->first] = matset;
545 
546  rval = mbImpl->add_entities( matset, mit->second );MB_CHK_SET_ERR( rval, "Couldn't add entities to material set" );
547  }
548 
549  return MB_SUCCESS;
550 }
551 
552 ErrorCode ReadCCMIO::construct_cells( TupleList& face_map,
553 #ifndef TUPLE_LIST
554  SenseList& sense_map,
555 #endif
556  TupleList& /* vert_map */,
557  std::map< int, int >& cell_topo_types,
558  std::vector< EntityHandle >& new_cells )
559 {
560  std::vector< EntityHandle > facehs;
561  std::vector< int > senses;
562  EntityHandle cell;
563  ErrorCode tmp_rval, rval = MB_SUCCESS;
564  EntityType this_type = MBMAXTYPE;
565  bool has_mid_nodes = false;
566 #ifdef TUPLE_LIST
567  unsigned int i = 0;
568  while( i < face_map.n )
569  {
570  // Pull out face handles bounding the same cell
571  facehs.clear();
572  int this_id = face_map.get_int( i );
573  unsigned int inext = i;
574  while( face_map.get_int( inext ) == this_id && inext <= face_map.n )
575  {
576  inext++;
577  EntityHandle face = face_map.get_ulong( inext );
578  facehs.push_back( face );
579  senses.push_back( face_map.get_short( inext ) );
580  }
581  this_type = MBMAXTYPE;
582  has_mid_nodes = false;
583 #else
584  std::map< int, std::vector< EntityHandle > >::iterator fmit;
585  std::map< int, std::vector< int > >::iterator smit;
586  std::map< int, int >::iterator typeit;
587  for( fmit = face_map.begin(), smit = sense_map.begin(); fmit != face_map.end(); ++fmit, ++smit )
588  {
589  // Pull out face handles bounding the same cell
590  facehs.clear();
591  int this_id = ( *fmit ).first;
592  facehs = ( *fmit ).second;
593  senses.clear();
594  senses = ( *smit ).second;
595  typeit = cell_topo_types.find( this_id );
596  if( typeit != cell_topo_types.end() )
597  {
598  rval = ccmio_to_moab_type( typeit->second, this_type, has_mid_nodes );
599  }
600  else
601  {
602  this_type = MBMAXTYPE;
603  has_mid_nodes = false;
604  }
605 #endif
606  tmp_rval = create_cell_from_faces( facehs, senses, this_type, has_mid_nodes, cell );
607  if( MB_SUCCESS != tmp_rval )
608  rval = tmp_rval;
609  else
610  {
611  new_cells.push_back( cell );
612  // Tag cell with global id
613  tmp_rval = mbImpl->tag_set_data( mGlobalIdTag, &cell, 1, &this_id );
614  if( MB_SUCCESS != tmp_rval ) rval = tmp_rval;
615  }
616  }
617 
618  return rval;
619 }
620 
621 ErrorCode ReadCCMIO::ccmio_to_moab_type( int ccm_type, EntityType& moab_type, bool& has_mid_nodes )
622 {
623  switch( ccm_type )
624  {
625  case 1:
626  moab_type = MBVERTEX;
627  break;
628  case 2:
629  case 28:
630  moab_type = MBEDGE;
631  break;
632  case 29:
633  moab_type = MBMAXTYPE;
634  break;
635  case 3:
636  case 4:
637  moab_type = MBQUAD;
638  break;
639  case 11:
640  case 21:
641  moab_type = MBHEX;
642  break;
643  case 12:
644  case 22:
645  moab_type = MBPRISM;
646  break;
647  case 13:
648  case 23:
649  moab_type = MBTET;
650  break;
651  case 14:
652  case 24:
653  moab_type = MBPYRAMID;
654  break;
655  case 255:
656  moab_type = MBPOLYHEDRON;
657  break;
658  default:
659  moab_type = MBMAXTYPE;
660  }
661 
662  switch( ccm_type )
663  {
664  case 28:
665  case 4:
666  case 21:
667  case 22:
668  case 23:
669  case 24:
670  has_mid_nodes = true;
671  break;
672  default:
673  break;
674  }
675 
676  return MB_SUCCESS;
677 }
678 
679 ErrorCode ReadCCMIO::create_cell_from_faces( std::vector< EntityHandle >& facehs,
680  std::vector< int >& senses,
681  EntityType this_type,
682  bool /* has_mid_nodes */,
683  EntityHandle& cell )
684 {
685  ErrorCode rval;
686 
687  // Test up front to see if they're one type
688  EntityType face_type = mbImpl->type_from_handle( facehs[0] );
689  bool same_type = true;
690  for( std::vector< EntityHandle >::iterator vit = facehs.begin(); vit != facehs.end(); ++vit )
691  {
692  if( face_type != mbImpl->type_from_handle( *vit ) )
693  {
694  same_type = false;
695  break;
696  }
697  }
698 
699  std::vector< EntityHandle > verts;
700  EntityType input_type = this_type;
701  std::vector< EntityHandle > storage;
702  MeshTopoUtil mtu( mbImpl );
703 
704  // Preset this to maxtype, so we get an affirmative choice in loop
705  this_type = MBMAXTYPE;
706 
707  if( ( MBTET == input_type || MBMAXTYPE == input_type ) && same_type && face_type == MBTRI && facehs.size() == 4 )
708  {
709  // Try to get proper connectivity for tet
710 
711  // Get connectivity of first face, and reverse it if sense is forward, since
712  // base face always points into entity
713  rval = mbImpl->get_connectivity( &facehs[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
714  if( senses[0] > 0 ) std::reverse( verts.begin(), verts.end() );
715 
716  // Get the 4th vertex through the next tri
717  const EntityHandle* conn;
718  int conn_size;
719  rval = mbImpl->get_connectivity( facehs[1], conn, conn_size, true, &storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
720  int i = 0;
721  while( std::find( verts.begin(), verts.end(), conn[i] ) != verts.end() && i < conn_size )
722  i++;
723 
724  // If i is not at the end of the verts, found the apex; otherwise fall back to polyhedron
725  if( conn_size != i )
726  {
727  this_type = MBTET;
728  verts.push_back( conn[i] );
729  }
730  }
731  else if( ( MBHEX == input_type || MBMAXTYPE == input_type ) && same_type && MBQUAD == face_type &&
732  facehs.size() == 6 )
733  {
734  // Build hex from quads
735  // Algorithm:
736  // - verts = vertices from 1st quad
737  // - Find quad q1 sharing verts[0] and verts[1]
738  // - Find quad q2 sharing other 2 verts in q1
739  // - Find v1 = opposite vert from verts[1] in q1 , v2 = opposite from verts[0]
740  // - Get i = offset of v1 in verts2 of q2, rotate verts2 by i
741  // - If verts2[(i + 1) % 4] != v2, flip verts2 by switching verts2[1] and verts2[3]
742  // - append verts2 to verts
743 
744  // Get the other vertices for this hex; need to find the quad with no common vertices
745  Range tmp_faces, tmp_verts;
746  // Get connectivity of first face, and reverse it if sense is forward, since
747  // base face always points into entity
748  rval = mbImpl->get_connectivity( &facehs[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
749  if( senses[0] > 0 ) std::reverse( verts.begin(), verts.end() );
750 
751  // Get q1, which shares 2 vertices with q0
752  std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
753  rval = mbImpl->get_adjacencies( &verts[0], 2, 2, false, tmp_faces );
754  if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
755  tmp_faces.erase( facehs[0] );
756  EntityHandle q1 = *tmp_faces.begin();
757  // Get other 2 verts of q1
758  rval = mbImpl->get_connectivity( &q1, 1, tmp_verts );MB_CHK_SET_ERR( rval, "Couldn't get adj verts" );
759  tmp_verts.erase( verts[0] );
760  tmp_verts.erase( verts[1] );
761  // Get q2
762  std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
763  rval = mbImpl->get_adjacencies( tmp_verts, 2, false, tmp_faces );
764  if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
765  tmp_faces.erase( q1 );
766  EntityHandle q2 = *tmp_faces.begin();
767  // Get verts in q2
768  rval = mbImpl->get_connectivity( &q2, 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get adj vertices" );
769 
770  // Get verts in q1 opposite from v[1] and v[0] in q0
771  EntityHandle v0 = 0, v1 = 0;
772  rval = mtu.opposite_entity( q1, verts[1], v0 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
773  rval = mtu.opposite_entity( q1, verts[0], v1 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
774  if( v0 && v1 )
775  {
776  // Offset of v0 in q2, then rotate and flip
777  unsigned int ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
778  if( 4 == ioff ) MB_SET_ERR( MB_FAILURE, "Trouble finding offset" );
779 
780  if( storage[( ioff + 1 ) % 4] != v1 )
781  {
782  std::reverse( storage.begin(), storage.end() );
783  ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
784  }
785  if( 0 != ioff ) std::rotate( storage.begin(), storage.begin() + ioff, storage.end() );
786 
787  // Copy into verts, and make hex
788  std::copy( storage.begin(), storage.end(), std::back_inserter( verts ) );
789  this_type = MBHEX;
790  }
791  }
792 
793  if( MBMAXTYPE == this_type && facehs.size() == 5 )
794  {
795  // Some preliminaries
796  std::vector< EntityHandle > tris, quads;
797  for( unsigned int i = 0; i < 5; i++ )
798  {
799  if( MBTRI == mbImpl->type_from_handle( facehs[i] ) )
800  tris.push_back( facehs[i] );
801  else if( MBQUAD == mbImpl->type_from_handle( facehs[i] ) )
802  quads.push_back( facehs[i] );
803  }
804 
805  // Check for prisms
806  if( 2 == tris.size() && 3 == quads.size() )
807  {
808  // OK, we have the right number of tris and quads; try to find the proper verts
809 
810  // Get connectivity of first tri, and reverse if necessary
811  int index = std::find( facehs.begin(), facehs.end(), tris[0] ) - facehs.begin();
812  rval = mbImpl->get_connectivity( &tris[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
813  if( senses[index] > 0 ) std::reverse( verts.begin(), verts.end() );
814 
815  // Now align vertices of other tri, through a quad, similar to how we did hexes
816  // Get q1, which shares 2 vertices with t0
817  Range tmp_faces, tmp_verts;
818  std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
819  rval = mbImpl->get_adjacencies( &verts[0], 2, 2, false, tmp_faces );
820  if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
821  tmp_faces.erase( tris[0] );
822  EntityHandle q1 = *tmp_faces.begin();
823  // Get verts in q1
824  rval = mbImpl->get_connectivity( &q1, 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get adj vertices" );
825 
826  // Get verts in q1 opposite from v[1] and v[0] in q0
827  EntityHandle v0 = 0, v1 = 0;
828  rval = mtu.opposite_entity( q1, verts[1], v0 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
829  rval = mtu.opposite_entity( q1, verts[0], v1 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
830  if( v0 && v1 )
831  {
832  // Offset of v0 in t2, then rotate and flip
833  storage.clear();
834  rval = mbImpl->get_connectivity( &tris[1], 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
835 
836  index = std::find( facehs.begin(), facehs.end(), tris[1] ) - facehs.begin();
837  if( senses[index] < 0 ) std::reverse( storage.begin(), storage.end() );
838  unsigned int ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
839  if( 3 == ioff ) MB_SET_ERR( MB_FAILURE, "Trouble finding offset" );
840  for( unsigned int i = 0; i < 3; i++ )
841  verts.push_back( storage[( ioff + i ) % 3] );
842 
843  this_type = MBPRISM;
844  }
845  }
846  else if( tris.size() == 4 && quads.size() == 1 )
847  {
848  // Check for pyramid
849  // Get connectivity of first tri, and reverse if necessary
850  int index = std::find( facehs.begin(), facehs.end(), quads[0] ) - facehs.begin();
851  rval = mbImpl->get_connectivity( &quads[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
852  if( senses[index] > 0 ) std::reverse( verts.begin(), verts.end() );
853 
854  // Get apex node
855  rval = mbImpl->get_connectivity( &tris[0], 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
856  for( unsigned int i = 0; i < 3; i++ )
857  {
858  if( std::find( verts.begin(), verts.end(), storage[i] ) == verts.end() )
859  {
860  verts.push_back( storage[i] );
861  break;
862  }
863  }
864 
865  if( 5 == verts.size() ) this_type = MBPYRAMID;
866  }
867  else
868  {
869  // Dummy else clause to stop in debugger
870  this_type = MBMAXTYPE;
871  }
872  }
873 
874  if( MBMAXTYPE != input_type && input_type != this_type && this_type != MBMAXTYPE )
875  std::cerr << "Warning: types disagree (cell_topo_type = " << CN::EntityTypeName( input_type )
876  << ", faces indicate type " << CN::EntityTypeName( this_type ) << std::endl;
877 
878  if( MBMAXTYPE != input_type && this_type == MBMAXTYPE && input_type != MBPOLYHEDRON )
879  std::cerr << "Warning: couldn't find proper connectivity for specified topo_type = "
880  << CN::EntityTypeName( input_type ) << std::endl;
881 
882  // Now make the element; if we fell back to polyhedron, use faces, otherwise use verts
883  if( MBPOLYHEDRON == this_type || MBMAXTYPE == this_type )
884  {
885  rval = mbImpl->create_element( MBPOLYHEDRON, &facehs[0], facehs.size(), cell );MB_CHK_SET_ERR( rval, "create_element failed" );
886  }
887  else
888  {
889  rval = mbImpl->create_element( this_type, &verts[0], verts.size(), cell );MB_CHK_SET_ERR( rval, "create_element failed" );
890  }
891 
892  return MB_SUCCESS;
893 }
894 
895 ErrorCode ReadCCMIO::read_all_faces( CCMIOID topologyID,
896  TupleList& vert_map,
897  TupleList& face_map,
898 #ifndef TUPLE_LIST
899  SenseList& sense_map,
900 #endif
901  Range* new_faces )
902 {
903  CCMIOSize_t index;
904  CCMIOID faceID;
905  ErrorCode rval;
906  CCMIOError error = kCCMIONoErr;
907 
908  // Get total # internal/bdy faces, size the face map accordingly
909 #ifdef TUPLE_LIST
910  index = CCMIOSIZEC( 0 );
911  int nbdy_faces = 0;
912  CCMIOSize_t nf;
913  error = kCCMIONoErr;
914  while( kCCMIONoErr == CCMIONextEntity( NULL, topologyID, kCCMIOBoundaryFaces, &index, &faceID ) )
915  {
916  CCMIOEntitySize( &error, faceID, &nf, NULL );
917  nbdy_faces += nf;
918  }
919  CCMIOGetEntity( &error, topologyID, kCCMIOInternalFaces, 0, &faceID );
920  CCMIOEntitySize( &error, faceID, &nf, NULL );
921 
922  int nint_faces = nf;
923  face_map.resize( 2 * nint_faces + nbdy_faces );
924 #endif
925 
926  // Get multiple blocks of bdy faces
927  index = CCMIOSIZEC( 0 );
928  while( kCCMIONoErr == CCMIONextEntity( NULL, topologyID, kCCMIOBoundaryFaces, &index, &faceID ) )
929  {
930  rval = read_faces( faceID, kCCMIOBoundaryFaces, vert_map, face_map,
931 #ifndef TUPLE_LIST
932  sense_map,
933 #endif
934  new_faces );MB_CHK_SET_ERR( rval, "Trouble reading boundary faces" );
935  }
936 
937  // Now get internal faces
938  CCMIOGetEntity( &error, topologyID, kCCMIOInternalFaces, 0, &faceID );
939  CHK_SET_CCMERR( error, "Couldn't get internal faces" );
940 
941  rval = read_faces( faceID, kCCMIOInternalFaces, vert_map, face_map,
942 #ifndef TUPLE_LIST
943  sense_map,
944 #endif
945  new_faces );MB_CHK_SET_ERR( rval, "Trouble reading internal faces" );
946 
947  return rval;
948 }
949 
950 ErrorCode ReadCCMIO::read_faces( CCMIOID faceID,
951  CCMIOEntity bdy_or_int,
952  TupleList& vert_map,
953  TupleList& face_map,
954 #ifndef TUPLE_LIST
955  SenseList& sense_map,
956 #endif
957  Range* new_faces )
958 {
959  if( kCCMIOInternalFaces != bdy_or_int && kCCMIOBoundaryFaces != bdy_or_int )
960  MB_SET_ERR( MB_FAILURE, "Face type isn't boundary or internal" );
961 
962  CCMIOSize_t dum_faces;
963  CCMIOError error = kCCMIONoErr;
964  CCMIOEntitySize( &error, faceID, &dum_faces, NULL );
965  int num_faces = GETINT32( dum_faces );
966 
967  // Get the size of the face connectivity array (not really a straight connect
968  // array, has n, connect(n), ...)
969  CCMIOSize_t farray_size = CCMIOSIZEC( 0 );
970  CCMIOReadFaces( &error, faceID, bdy_or_int, NULL, &farray_size, NULL, CCMIOINDEXC( kCCMIOStart ),
971  CCMIOINDEXC( kCCMIOEnd ) );
972  CHK_SET_CCMERR( error, "Trouble reading face connectivity length" );
973 
974  // Allocate vectors for holding farray and cells for each face; use new for finer
975  // control of de-allocation
976  int num_sides = ( kCCMIOInternalFaces == bdy_or_int ? 2 : 1 );
977  int* farray = new int[GETINT32( farray_size )];
978 
979  // Read farray and make the faces
980  CCMIOID mapID;
981  CCMIOReadFaces( &error, faceID, bdy_or_int, &mapID, NULL, farray, CCMIOINDEXC( kCCMIOStart ),
982  CCMIOINDEXC( kCCMIOEnd ) );
983  CHK_SET_CCMERR( error, "Trouble reading face connectivity" );
984 
985  std::vector< EntityHandle > face_handles;
986  ErrorCode rval = make_faces( farray, vert_map, face_handles, num_faces );MB_CHK_SET_ERR( rval, "Failed to make the faces" );
987 
988  // Read face cells and make tuples
989  int* face_cells;
990  if( num_sides * num_faces < farray_size )
991  face_cells = new int[num_sides * num_faces];
992  else
993  face_cells = farray;
994  CCMIOReadFaceCells( &error, faceID, bdy_or_int, face_cells, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
995  CHK_SET_CCMERR( error, "Trouble reading face cells" );
996 
997  int* tmp_ptr = face_cells;
998  for( unsigned int i = 0; i < face_handles.size(); i++ )
999  {
1000 #ifdef TUPLE_LIST
1001  short forward = 1, reverse = -1;
1002  face_map.push_back( &forward, tmp_ptr++, &face_handles[i], NULL );
1003  if( 2 == num_sides ) face_map.push_back( &reverse, tmp_ptr++, &face_handles[i], NULL );
1004 #else
1005  face_map[*tmp_ptr].push_back( face_handles[i] );
1006  sense_map[*tmp_ptr++].push_back( 1 );
1007  if( 2 == num_sides )
1008  {
1009  face_map[*tmp_ptr].push_back( face_handles[i] );
1010  sense_map[*tmp_ptr++].push_back( -1 );
1011  }
1012 #endif
1013  }
1014 
1015  // Now read & set face gids, reuse face_cells 'cuz we know it's big enough
1016  CCMIOReadMap( &error, mapID, face_cells, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
1017  CHK_SET_CCMERR( error, "Trouble reading face gids" );
1018 
1019  rval = mbImpl->tag_set_data( mGlobalIdTag, &face_handles[0], face_handles.size(), face_cells );MB_CHK_SET_ERR( rval, "Couldn't set face global ids" );
1020 
1021  // Make a neumann set for these faces if they're all in a boundary face set
1022  if( kCCMIOBoundaryFaces == bdy_or_int )
1023  {
1024  EntityHandle neuset;
1025  rval = mbImpl->create_meshset( MESHSET_SET, neuset );MB_CHK_SET_ERR( rval, "Failed to create neumann set" );
1026 
1027  // Don't trust entity index passed in
1028  int index;
1029  CCMIOGetEntityIndex( &error, faceID, &index );
1030  newNeusets[index] = neuset;
1031 
1032  rval = mbImpl->add_entities( neuset, &face_handles[0], face_handles.size() );MB_CHK_SET_ERR( rval, "Failed to add faces to neumann set" );
1033 
1034  // Now tag as neumann set; will add id later
1035  int dum_val = 0;
1036  rval = mbImpl->tag_set_data( mNeumannSetTag, &neuset, 1, &dum_val );MB_CHK_SET_ERR( rval, "Failed to tag neumann set" );
1037  }
1038 
1039  if( new_faces )
1040  {
1041  std::sort( face_handles.begin(), face_handles.end() );
1042  std::copy( face_handles.rbegin(), face_handles.rend(), range_inserter( *new_faces ) );
1043  }
1044 
1045  return MB_SUCCESS;
1046 }
1047 
1048 ErrorCode ReadCCMIO::make_faces( int* farray,
1049  TupleList& vert_map,
1050  std::vector< EntityHandle >& new_faces,
1051  int num_faces )
1052 {
1053  std::vector< EntityHandle > verts;
1054  ErrorCode tmp_rval = MB_SUCCESS, rval = MB_SUCCESS;
1055 
1056  for( int i = 0; i < num_faces; i++ )
1057  {
1058  int num_verts = *farray++;
1059  verts.resize( num_verts );
1060 
1061  // Fill in connectivity by looking up by gid in vert tuple_list
1062  for( int j = 0; j < num_verts; j++ )
1063  {
1064 #ifdef TUPLE_LIST
1065  int tindex = vert_map.find( 1, farray[j] );
1066  if( -1 == tindex )
1067  {
1068  tmp_rval = MB_FAILURE;
1069  break;
1070  }
1071  verts[j] = vert_map.get_ulong( tindex, 0 );
1072 #else
1073  verts[j] = ( vert_map[farray[j]] )[0];
1074 #endif
1075  }
1076  farray += num_verts;
1077 
1078  if( MB_SUCCESS == tmp_rval )
1079  {
1080  // Make face
1081  EntityType ftype = ( 3 == num_verts ? MBTRI : ( 4 == num_verts ? MBQUAD : MBPOLYGON ) );
1082  EntityHandle faceh;
1083  tmp_rval = mbImpl->create_element( ftype, &verts[0], num_verts, faceh );
1084  if( faceh ) new_faces.push_back( faceh );
1085  }
1086 
1087  if( MB_SUCCESS != tmp_rval ) rval = tmp_rval;
1088  }
1089 
1090  return rval;
1091 }
1092 
1093 ErrorCode ReadCCMIO::read_vertices( CCMIOSize_t /* proc */,
1094  CCMIOID /* processorID */,
1095  CCMIOID verticesID,
1096  CCMIOID /* topologyID */,
1097  Range* verts,
1098  TupleList& vert_map )
1099 {
1100  CCMIOError error = kCCMIONoErr;
1101 
1102  // Pre-read the number of vertices, so we can pre-allocate & read directly in
1103  CCMIOSize_t nverts = CCMIOSIZEC( 0 );
1104  CCMIOEntitySize( &error, verticesID, &nverts, NULL );
1105  CHK_SET_CCMERR( error, "Couldn't get number of vertices" );
1106 
1107  // Get # dimensions
1108  CCMIOSize_t dims;
1109  float scale;
1110  CCMIOReadVerticesf( &error, verticesID, &dims, NULL, NULL, NULL, CCMIOINDEXC( 0 ), CCMIOINDEXC( 1 ) );
1111  CHK_SET_CCMERR( error, "Couldn't get number of dimensions" );
1112 
1113  // Allocate vertex space
1114  EntityHandle node_handle = 0;
1115  std::vector< double* > arrays;
1116  readMeshIface->get_node_coords( 3, GETINT32( nverts ), MB_START_ID, node_handle, arrays );
1117 
1118  // Read vertex coords
1119  CCMIOID mapID;
1120  std::vector< double > tmp_coords( GETINT32( dims ) * GETINT32( nverts ) );
1121  CCMIOReadVerticesd( &error, verticesID, &dims, &scale, &mapID, &tmp_coords[0], CCMIOINDEXC( 0 ),
1122  CCMIOINDEXC( 0 + nverts ) );
1123  CHK_SET_CCMERR( error, "Trouble reading vertex coordinates" );
1124 
1125  // Copy interleaved coords into moab blocked coordinate space
1126  int i = 0, threei = 0;
1127  for( ; i < nverts; i++ )
1128  {
1129  arrays[0][i] = tmp_coords[threei++];
1130  arrays[1][i] = tmp_coords[threei++];
1131  if( 3 == GETINT32( dims ) )
1132  arrays[2][i] = tmp_coords[threei++];
1133  else
1134  arrays[2][i] = 0.0;
1135  }
1136 
1137  // Scale, if necessary
1138  if( 1.0 != scale )
1139  {
1140  for( i = 0; i < nverts; i++ )
1141  {
1142  arrays[0][i] *= scale;
1143  arrays[1][i] *= scale;
1144  if( 3 == GETINT32( dims ) ) arrays[2][i] *= scale;
1145  }
1146  }
1147 
1148  // Read gids for vertices
1149  std::vector< int > gids( GETINT32( nverts ) );
1150  CCMIOReadMap( &error, mapID, &gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
1151  CHK_SET_CCMERR( error, "Trouble reading vertex global ids" );
1152 
1153  // Put new vertex handles into range, and set gids for them
1154  Range new_verts( node_handle, node_handle + nverts - 1 );
1155  ErrorCode rval = mbImpl->tag_set_data( mGlobalIdTag, new_verts, &gids[0] );MB_CHK_SET_ERR( rval, "Couldn't set gids on vertices" );
1156 
1157  // Pack vert_map with global ids and handles for these vertices
1158 #ifdef TUPLE_LIST
1159  vert_map.resize( GETINT32( nverts ) );
1160  for( i = 0; i < GETINT32( nverts ); i++ )
1161  {
1162  vert_map.push_back( NULL, &gids[i], &node_handle, NULL );
1163 #else
1164  for( i = 0; i < GETINT32( nverts ); i++ )
1165  {
1166  ( vert_map[gids[i]] ).push_back( node_handle );
1167 #endif
1168  node_handle += 1;
1169  }
1170 
1171  if( verts ) verts->merge( new_verts );
1172 
1173  return MB_SUCCESS;
1174 }
1175 
1176 ErrorCode ReadCCMIO::get_processors( CCMIOID stateID,
1177  CCMIOID& processorID,
1178  CCMIOID& verticesID,
1179  CCMIOID& topologyID,
1180  CCMIOID& solutionID,
1181  std::vector< CCMIOSize_t >& procs,
1182  bool& /* has_solution */ )
1183 {
1184  CCMIOSize_t proc = CCMIOSIZEC( 0 );
1185  CCMIOError error = kCCMIONoErr;
1186 
1187  CCMIONextEntity( &error, stateID, kCCMIOProcessor, &proc, &processorID );
1188  CHK_SET_CCMERR( error, "CCMIONextEntity() failed" );
1189  if( CCMIOReadProcessor( NULL, processorID, &verticesID, &topologyID, NULL, &solutionID ) != kCCMIONoErr )
1190  {
1191  // Maybe no solution; try again
1192  CCMIOReadProcessor( &error, processorID, &verticesID, &topologyID, NULL, NULL );
1193  CHK_SET_CCMERR( error, "CCMIOReadProcessor() failed" );
1194  hasSolution = false;
1195  }
1196 
1197  procs.push_back( proc );
1198 
1199  return MB_SUCCESS;
1200 }
1201 
1202 ErrorCode ReadCCMIO::read_tag_values( const char* /* file_name */,
1203  const char* /* tag_name */,
1204  const FileOptions& /* opts */,
1205  std::vector< int >& /* tag_values_out */,
1206  const SubsetList* /* subset_list */ )
1207 {
1208  return MB_FAILURE;
1209 }
1210 
1211 } // namespace moab