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