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