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