1 /** \file MsqMOAB.cpp 2 * \brief 3 * \author Vijay Mahadevan 4 */ 5 #include "MsqMOAB.hpp" 6 #include "MsqError.hpp" 7 #include "MeshInterface.hpp" 8 #include "MsqDebug.hpp" 9 #include "MsqVertex.hpp" 10 #include <assert.h> 11 #include "MsqIBase.hpp" 12 #include <algorithm> 13 14 #ifdef IMESH_MAJOR_VERSION 15 # define IMESH_VERSION_ATLEAST(MAJOR,MINOR) \ 16 1000*IMESH_MAJOR_VERSION+IMESH_MINOR_VERSION <= \ 17 1000*MAJOR+MINOR 18 #else 19 # define IMESH_VERSION_ATLEAST(MAJOR,MINOR) 0 20 #endif 21 22 namespace MBMesquite 23 { 24 25 26 /************************************************************************* 27 * Mesh Definition 28 ************************************************************************/ 29 MsqMOAB(moab::Core * mesh,moab::EntityHandle meshset,moab::EntityType type,MsqError & err,const moab::Tag * fixed_tag,const moab::Tag * slaved_tag)30 MsqMOAB::MsqMOAB ( moab::Core* mesh, 31 moab::EntityHandle meshset, 32 moab::EntityType type, 33 MsqError& err, 34 const moab::Tag* fixed_tag, 35 const moab::Tag* slaved_tag ) 36 : meshInstance ( mesh ), 37 inputSetType ( moab::MBMAXTYPE ), 38 inputSet ( 0 ), 39 byteTag ( 0 ), 40 createdByteTag ( false ), 41 geometricDimension ( 0 ) 42 { 43 init_active_mesh ( mesh, err, fixed_tag, slaved_tag ); 44 MSQ_ERRRTN ( err ); 45 set_active_set ( meshset, type, err ); 46 MSQ_ERRRTN ( err ); 47 } 48 MsqMOAB(moab::Core * mesh,moab::EntityType type,MsqError & err,const moab::Tag * fixed_tag,const moab::Tag * slaved_tag)49 MsqMOAB::MsqMOAB ( moab::Core* mesh, 50 moab::EntityType type, 51 MsqError& err, 52 const moab::Tag* fixed_tag, 53 const moab::Tag* slaved_tag ) 54 : meshInstance ( mesh ), 55 inputSetType ( moab::MBMAXTYPE ), 56 inputSet ( 0 ), 57 byteTag ( 0 ), 58 createdByteTag ( false ), 59 geometricDimension ( 0 ) 60 { 61 init_active_mesh ( mesh, err, fixed_tag, slaved_tag ); 62 MSQ_ERRRTN ( err ); 63 64 moab::EntityHandle root_set = 0 /*meshInstance->get_root_set()*/; 65 set_active_set ( root_set, type, err ); MSQ_ERRRTN ( err ); 66 } 67 ~MsqMOAB()68 MsqMOAB::~MsqMOAB() 69 { 70 moab::ErrorCode ierr; 71 72 if ( createdByteTag ) 73 { 74 ierr = meshInstance->tag_delete ( byteTag ); MB_CHK_ERR_RET ( ierr ); 75 } 76 } 77 get_interface() const78 moab::Core* MsqMOAB::get_interface() const 79 { 80 return meshInstance; 81 } 82 get_entity_set() const83 moab::EntityHandle MsqMOAB::get_entity_set() const 84 { 85 return inputSet; 86 } 87 check_valid_flag_tag(moab::Tag tag,const char *,MsqError & err)88 moab::DataType MsqMOAB::check_valid_flag_tag ( moab::Tag tag, 89 const char* /*which_flag*/, 90 MsqError& err ) 91 { 92 moab::ErrorCode ierr; 93 int size; 94 std::string name; 95 moab::DataType type = moab::MB_MAX_DATA_TYPE; 96 97 ierr = meshInstance->tag_get_data_type ( tag, type ); MB_CHK_ERR_RET_VAL ( ierr, type ); 98 ierr = meshInstance->tag_get_name ( tag, name ); MB_CHK_ERR_RET_VAL ( ierr, type ); 99 ierr = meshInstance->tag_get_length ( tag, size ); MB_CHK_ERR_RET_VAL ( ierr, type ); 100 101 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 102 return type; 103 } 104 init_active_mesh(moab::Core *,MsqError & err,const moab::Tag * fixed_tag,const moab::Tag * slaved_tag)105 void MsqMOAB::init_active_mesh ( moab::Core* /*mesh*/, 106 MsqError& err, 107 const moab::Tag* fixed_tag, 108 const moab::Tag* slaved_tag ) 109 { 110 moab::ErrorCode ierr; 111 112 // Initialize topology map 113 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 114 115 const size_t mapsize = sizeof ( topologyMap ) / sizeof ( MBMesquite::EntityTopology ); 116 117 if ( mapsize < moab::MBMAXTYPE ) 118 { 119 MSQ_SETERR ( err ) ( "MsqMOAB needs to be updated for new iMesh element topologies.", 120 MsqError::INTERNAL_ERROR ); 121 } 122 123 for ( size_t i = 0; i <= moab::MBMAXTYPE; ++i ) 124 { topologyMap[i] = MBMesquite::MIXED; } 125 126 topologyMap[moab::MBTRI ] = MBMesquite::TRIANGLE; 127 topologyMap[moab::MBQUAD ] = MBMesquite::QUADRILATERAL; 128 topologyMap[moab::MBTET ] = MBMesquite::TETRAHEDRON; 129 topologyMap[moab::MBHEX ] = MBMesquite::HEXAHEDRON; 130 topologyMap[moab::MBPRISM ] = MBMesquite::PRISM; 131 topologyMap[moab::MBPYRAMID ] = MBMesquite::PYRAMID; 132 133 // Check that fixed tag is valid 134 haveFixedTag = false; 135 136 if ( fixed_tag ) 137 { 138 fixedTagType = check_valid_flag_tag ( *fixed_tag, "fixed", err ); 139 MSQ_ERRRTN ( err ); 140 haveFixedTag = true; 141 fixedTag = *fixed_tag; 142 } 143 144 // Check that slaved tag is valid 145 haveSlavedTag = false; 146 147 if ( slaved_tag ) 148 { 149 slavedTagType = check_valid_flag_tag ( *slaved_tag, "slaved", err ); 150 MSQ_ERRRTN ( err ); 151 haveSlavedTag = true; 152 slavedTag = *slaved_tag; 153 } 154 155 // Get/create tag for vertex byte 156 ierr = meshInstance->tag_get_handle ( VERTEX_BYTE_TAG_NAME, byteTag ); 157 158 if ( moab::MB_SUCCESS != ierr ) 159 { 160 int defval = 0; 161 ierr = meshInstance->tag_get_handle ( VERTEX_BYTE_TAG_NAME, 1, moab::MB_TYPE_INTEGER, byteTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &defval, &createdByteTag ); //MB_CHK_ERR_RET(ierr); 162 } 163 else 164 { 165 int size; 166 moab::DataType type; 167 ierr = meshInstance->tag_get_length ( byteTag, size ); MB_CHK_ERR_RET ( ierr ); 168 ierr = meshInstance->tag_get_data_type ( byteTag, type ); MB_CHK_ERR_RET ( ierr ); 169 } 170 171 ierr = meshInstance->get_dimension ( geometricDimension ); MB_CHK_ERR_RET ( ierr ); 172 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 173 } 174 175 set_fixed_tag(moab::Tag tag,MsqError & err)176 void MsqMOAB::set_fixed_tag ( moab::Tag tag, MsqError& err ) 177 { 178 moab::DataType t = check_valid_flag_tag ( tag, "fixed", err ); 179 MSQ_ERRRTN ( err ); 180 fixedTag = tag; 181 fixedTagType = t; 182 haveFixedTag = true; 183 } 184 clear_fixed_tag()185 void MsqMOAB::clear_fixed_tag() 186 { 187 haveFixedTag = false; 188 } 189 get_fixed_tag() const190 const moab::Tag* MsqMOAB::get_fixed_tag() const 191 { 192 return haveFixedTag ? &fixedTag : 0; 193 } 194 set_slaved_tag(moab::Tag tag,MsqError & err)195 void MsqMOAB::set_slaved_tag ( moab::Tag tag, MsqError& err ) 196 { 197 moab::DataType t = check_valid_flag_tag ( tag, "slaved", err ); 198 MSQ_ERRRTN ( err ); 199 slavedTag = tag; 200 slavedTagType = t; 201 haveSlavedTag = true; 202 } 203 clear_slaved_tag()204 void MsqMOAB::clear_slaved_tag() 205 { 206 haveSlavedTag = false; 207 } 208 get_slaved_tag() const209 const moab::Tag* MsqMOAB::get_slaved_tag() const 210 { 211 return haveSlavedTag ? &slavedTag : 0; 212 } 213 214 215 set_active_set(moab::EntityHandle elem_set,moab::EntityType type_in,MsqError & err)216 void MsqMOAB::set_active_set ( moab::EntityHandle elem_set, 217 moab::EntityType type_in, 218 MsqError& err ) 219 { 220 inputSetType = type_in; 221 inputSet = elem_set; 222 223 // clear vertex byte 224 std::vector<VertexHandle> verts; 225 get_all_vertices ( verts, err ); MSQ_ERRRTN ( err ); 226 227 if ( !verts.empty() ) 228 { 229 std::vector<unsigned char> zeros ( verts.size(), 0 ); 230 vertices_set_byte ( arrptr ( verts ), arrptr ( zeros ), verts.size(), err ); 231 MSQ_CHKERR ( err ); 232 } 233 } 234 235 236 237 // Returns whether this mesh lies in a 2D or 3D coordinate system. get_geometric_dimension(MBMesquite::MsqError &)238 int MsqMOAB::get_geometric_dimension ( MBMesquite::MsqError& /*err*/ ) 239 { 240 return geometricDimension; 241 } 242 243 244 //************ Vertex Properties ******************** 245 get_flag_data(moab::Tag tag,bool have_tag,moab::DataType type,const VertexHandle vert_array[],std::vector<bool> & flag_array,size_t num_vtx,MsqError & err)246 void MsqMOAB::get_flag_data ( moab::Tag tag, 247 bool have_tag, 248 moab::DataType type, 249 const VertexHandle vert_array[], 250 std::vector<bool>& flag_array, 251 size_t num_vtx, 252 MsqError& err ) 253 { 254 if ( !num_vtx ) 255 { return; } 256 257 if ( !have_tag ) 258 { 259 flag_array.clear(); 260 flag_array.resize ( num_vtx, false ); 261 return; 262 } 263 264 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 265 flag_array.resize ( num_vtx ); 266 267 assert ( sizeof ( VertexHandle ) == sizeof ( moab::EntityHandle ) ); 268 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( vert_array ); 269 270 moab::ErrorCode ierr; 271 int alloc = num_vtx; 272 assert ( ( size_t ) alloc == num_vtx ); // size_t can hold larger values than int if 64-bit 273 274 if ( type == moab::MB_TYPE_INTEGER ) 275 { 276 std::vector<int> values ( num_vtx ); 277 int* ptr = arrptr ( values ); 278 ierr = meshInstance->tag_get_data ( tag, arr, num_vtx, ptr ); MB_CHK_ERR_RET ( ierr ); 279 280 for ( size_t i = 0; i < num_vtx; ++i ) 281 { flag_array[i] = !!values[i]; } 282 } 283 else if ( type == moab::MB_TYPE_OPAQUE ) 284 { 285 std::vector<char> values ( num_vtx ); 286 void* ptr = arrptr ( values ); 287 ierr = meshInstance->tag_get_data ( tag, arr, num_vtx, ptr ); MB_CHK_ERR_RET ( ierr ); 288 289 for ( size_t i = 0; i < num_vtx; ++i ) 290 { flag_array[i] = !!values[i]; } 291 } 292 else 293 { 294 MSQ_SETERR ( err ) ( "Invalid tag type for vertex flag data", MsqError::INVALID_STATE ); 295 return ; 296 } 297 298 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 299 } 300 301 // Returns true or false, indicating whether the vertex 302 // is allowed to be repositioned. True indicates that the vertex 303 // is fixed and cannot be moved. Note that this is a read-only 304 // property; this flag can't be modified by users of the 305 // MBMesquite::Mesh interface. vertices_get_fixed_flag(const VertexHandle vert_array[],std::vector<bool> & bool_array,size_t num_vtx,MsqError & err)306 void MsqMOAB::vertices_get_fixed_flag ( 307 const VertexHandle vert_array[], 308 std::vector<bool>& bool_array, 309 size_t num_vtx, MsqError& err ) 310 { 311 get_flag_data ( fixedTag, haveFixedTag, fixedTagType, vert_array, bool_array, num_vtx, err ); 312 } 313 314 vertices_get_slaved_flag(const VertexHandle vert_array[],std::vector<bool> & bool_array,size_t num_vtx,MsqError & err)315 void MsqMOAB::vertices_get_slaved_flag ( 316 const VertexHandle vert_array[], 317 std::vector<bool>& bool_array, 318 size_t num_vtx, MsqError& err ) 319 { 320 get_flag_data ( slavedTag, haveSlavedTag, slavedTagType, vert_array, bool_array, num_vtx, err ); 321 } 322 323 // Get vertex coordinates vertices_get_coordinates(const MBMesquite::Mesh::VertexHandle vert_array[],MsqVertex * coordinates,size_t num_vtx,MsqError & err)324 void MsqMOAB::vertices_get_coordinates ( 325 const MBMesquite::Mesh::VertexHandle vert_array[], 326 MsqVertex* coordinates, 327 size_t num_vtx, 328 MsqError& err ) 329 { 330 if ( !num_vtx ) 331 { return; } 332 333 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 334 std::vector<double> dbl_store ( 3 * num_vtx ); 335 double* dbl_array = arrptr ( dbl_store ); 336 337 moab::ErrorCode ierr; 338 // int junk = 3*num_vtx; 339 assert ( sizeof ( VertexHandle ) == sizeof ( moab::EntityHandle ) ); 340 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( vert_array ); 341 ierr = meshInstance->get_coords ( arr, num_vtx, dbl_array ); MB_CHK_ERR_RET ( ierr ); 342 343 if ( geometricDimension == 2 ) 344 { 345 double* iter = dbl_array; 346 347 for ( size_t i = 0; i < num_vtx; ++i ) 348 { 349 coordinates[i].x ( *iter ); ++iter; 350 coordinates[i].y ( *iter ); ++iter; 351 coordinates[i].z ( 0 ); 352 } 353 } 354 else 355 { 356 double* iter = dbl_array; 357 358 for ( size_t i = 0; i < num_vtx; ++i ) 359 { 360 coordinates[i].set ( iter ); 361 iter += 3; 362 } 363 } 364 365 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 366 } 367 vertex_set_coordinates(MBMesquite::Mesh::VertexHandle vertex,const Vector3D & coords,MsqError & err)368 void MsqMOAB::vertex_set_coordinates ( 369 MBMesquite::Mesh::VertexHandle vertex, 370 const Vector3D& coords, MsqError& err ) 371 { 372 moab::ErrorCode ierr; 373 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 374 moab::EntityHandle* bh = reinterpret_cast<moab::EntityHandle*> ( &vertex ); 375 double xyz[3] = {coords[0], coords[1], coords[2]}; 376 ierr = meshInstance->set_coords ( bh, 1, xyz ); MB_CHK_ERR_RET ( ierr ); 377 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 378 } 379 380 // Each vertex has a byte-sized flag that can be used to store 381 // flags. This byte's value is neither set nor used by the mesh 382 // implementation. It is intended to be used by Mesquite algorithms. 383 // Until a vertex's byte has been explicitly set, its value is 0. vertex_set_byte(MBMesquite::Mesh::VertexHandle vertex,unsigned char byte,MsqError & err)384 void MsqMOAB::vertex_set_byte ( 385 MBMesquite::Mesh::VertexHandle vertex, 386 unsigned char byte, MsqError& err ) 387 { 388 moab::ErrorCode ierr; 389 int value = byte; 390 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 391 moab::EntityHandle* bh = reinterpret_cast<moab::EntityHandle*> ( &vertex ); 392 ierr = meshInstance->tag_set_data ( byteTag, bh, 1, &value ); MB_CHK_ERR_RET ( ierr ); 393 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 394 } 395 vertices_set_byte(const VertexHandle * vert_array,const unsigned char * byte_array,size_t array_size,MsqError & err)396 void MsqMOAB::vertices_set_byte ( 397 const VertexHandle* vert_array, 398 const unsigned char* byte_array, 399 size_t array_size, MsqError& err ) 400 { 401 if ( !array_size ) 402 { return; } 403 404 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 405 std::vector<int> data ( array_size ); 406 std::copy ( byte_array, byte_array + array_size, data.begin() ); 407 moab::ErrorCode ierr; 408 assert ( sizeof ( VertexHandle ) == sizeof ( moab::EntityHandle ) ); 409 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( vert_array ); 410 ierr = meshInstance->tag_set_data ( byteTag, arr, array_size, arrptr ( data ) ); MB_CHK_ERR_RET ( ierr ); 411 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 412 } 413 414 // Retrieve the byte value for the specified vertex or vertices. 415 // The byte value is 0 if it has not yet been set via one of the 416 // *_set_byte() functions. vertex_get_byte(MBMesquite::Mesh::VertexHandle vertex,unsigned char * byte,MsqError & err)417 void MsqMOAB::vertex_get_byte ( 418 MBMesquite::Mesh::VertexHandle vertex, 419 unsigned char* byte, MsqError& err ) 420 { 421 moab::ErrorCode ierr; 422 int value; 423 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 424 moab::EntityHandle* bh = reinterpret_cast<moab::EntityHandle*> ( &vertex ); 425 ierr = meshInstance->tag_get_data ( byteTag, bh, 1, &value ); MB_CHK_ERR_RET ( ierr ); 426 *byte = value; 427 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 428 } 429 vertices_get_byte(const VertexHandle * vert_array,unsigned char * byte_array,size_t array_size,MsqError & err)430 void MsqMOAB::vertices_get_byte ( 431 const VertexHandle* vert_array, 432 unsigned char* byte_array, 433 size_t array_size, MsqError& err ) 434 { 435 if ( !array_size ) 436 { return; } 437 438 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 439 std::vector<int> data ( array_size ); 440 moab::ErrorCode ierr; 441 int* ptr = arrptr ( data ); 442 assert ( sizeof ( VertexHandle ) == sizeof ( moab::EntityHandle ) ); 443 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( vert_array ); 444 ierr = meshInstance->tag_get_data ( byteTag, arr, array_size, ptr ); MB_CHK_ERR_RET ( ierr ); 445 std::copy ( data.begin(), data.end(), byte_array ); 446 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 447 } 448 449 450 //**************** Topology ***************** 451 get_adjacent_entities(const moab::EntityHandle * source,size_t num_source,moab::EntityType target_type,std::vector<EntityHandle> & target,std::vector<size_t> & offsets,MsqError & err)452 void MsqMOAB::get_adjacent_entities ( const moab::EntityHandle* source, 453 size_t num_source, 454 moab::EntityType target_type, 455 std::vector<EntityHandle>& target, 456 std::vector<size_t>& offsets, 457 MsqError& err ) 458 { 459 if ( num_source == 0 ) 460 { 461 target.clear(); 462 offsets.clear(); 463 offsets.reserve ( 1 ); 464 offsets.push_back ( 0 ); 465 return; 466 } 467 468 err.set_error ( MBMesquite::MsqError::UNKNOWN_ERROR ); 469 470 moab::ErrorCode ierr; 471 int num_adj = 0, num_offset = 0; 472 unsigned iadjoff = 0; 473 474 assert ( sizeof ( size_t ) >= sizeof ( int ) ); 475 offsets.resize ( num_source + 1 ); 476 int* ptr2; 477 bool expand = false; 478 479 if ( sizeof ( size_t ) > sizeof ( int ) ) 480 { 481 ptr2 = ( int* ) malloc ( sizeof ( int ) * ( num_source + 1 ) ); 482 expand = true; 483 num_offset = num_source + 1; 484 } 485 else 486 { 487 // sizeof(int) == sizeof(size_t) 488 ptr2 = ( int* ) arrptr ( offsets ); 489 num_offset = offsets.size(); 490 } 491 492 // std::cout << "Target capacity: " << target.capacity() << " and num sources = " << num_source << std::endl; 493 494 assert ( sizeof ( moab::EntityHandle ) == sizeof ( EntityHandle ) ); 495 bool have_adj = false; 496 497 // If passed vector has allocated storage, try to use existing space 498 if ( target.capacity() >= num_source ) 499 { 500 target.resize ( target.capacity() ); 501 // target.clear(); 502 moab::EntityHandle* ptr = reinterpret_cast<moab::EntityHandle*> ( arrptr ( target ) ); 503 504 ptr2[0] = 0; 505 506 for ( unsigned is = 0; is < num_source; ++is ) 507 { 508 moab::Range adjents; 509 510 if ( target_type == moab::MBVERTEX ) 511 { 512 ierr = meshInstance->get_adjacencies ( &source[is], 1, 0, false, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 513 } 514 else if ( target_type == moab::MBEDGE ) 515 { 516 ierr = meshInstance->get_adjacencies ( &source[is], 1, 1, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 517 } 518 else if ( target_type <= moab::MBPOLYGON ) 519 { 520 ierr = meshInstance->get_adjacencies ( &source[is], 1, 2, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 521 } 522 else if ( target_type < moab::MBENTITYSET ) 523 { 524 ierr = meshInstance->get_adjacencies ( &source[is], 1, 3, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 525 } 526 else /* Either EntitySet or MBMaxType -- Failures */ 527 { 528 MSQ_SETERR ( err ) ( process_itaps_error ( 1 ), MsqError::NOT_IMPLEMENTED ); // "Invalid Target entity type specified" 529 return; 530 } 531 532 ptr2[is + 1] = iadjoff + adjents.size(); 533 534 for ( unsigned iadj = 0; iadj < adjents.size(); ++iadj, ++iadjoff ) 535 { ptr[iadjoff] = adjents[iadj]; } 536 537 // target.push_back( static_cast<EntityHandle>(&adjents[iadj]) ); 538 // std::cout << "1. Source entity [ " << is << "]: n(adjacencies) = " << offsets[is+1] << std::endl; 539 } 540 541 if ( moab::MB_SUCCESS == ierr ) 542 { 543 have_adj = true; 544 num_adj = iadjoff; 545 } 546 } 547 548 // If implementation passed back a size, try that 549 if ( !have_adj && num_adj && ( unsigned ) num_adj > target.capacity() ) 550 { 551 target.resize ( num_adj ); 552 // int junk1 = target.capacity(), junk3 = offsets.size(); 553 moab::EntityHandle* ptr = reinterpret_cast<moab::EntityHandle*> ( arrptr ( target ) ); 554 555 ptr2[0] = 0; 556 557 for ( unsigned is = 0; is < num_source; ++is ) 558 { 559 moab::Range adjents; 560 561 if ( target_type == moab::MBVERTEX ) 562 { 563 ierr = meshInstance->get_adjacencies ( &source[is], 1, 0, false, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 564 } 565 else if ( target_type == moab::MBEDGE ) 566 { 567 ierr = meshInstance->get_adjacencies ( &source[is], 1, 1, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 568 } 569 else if ( target_type <= moab::MBPOLYGON ) 570 { 571 ierr = meshInstance->get_adjacencies ( &source[is], 1, 2, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 572 } 573 else if ( target_type < moab::MBENTITYSET ) 574 { 575 ierr = meshInstance->get_adjacencies ( &source[is], 1, 3, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 576 } 577 else /* Either EntitySet or MBMaxType -- Failures */ 578 { 579 MSQ_SETERR ( err ) ( process_itaps_error ( 1 ), MsqError::NOT_IMPLEMENTED ); // "Invalid Target entity type specified" 580 return; 581 } 582 583 ptr2[is + 1] = iadjoff + adjents.size(); 584 585 for ( unsigned iadj = 0; iadj < adjents.size(); ++iadj, ++iadjoff ) 586 { ptr[iadjoff] = adjents[iadj]; } 587 } 588 589 if ( moab::MB_SUCCESS == ierr ) 590 { 591 have_adj = true; 592 } 593 } 594 595 // Try with empty sidl array, and copy into elements vector 596 if ( !have_adj ) 597 { 598 // If implementation passed back a size, try that 599 std::vector<moab::EntityHandle> mArray; 600 601 ptr2[0] = iadjoff; 602 603 for ( unsigned is = 0; is < num_source; ++is ) 604 { 605 moab::Range adjents; 606 607 if ( target_type == moab::MBVERTEX ) 608 { 609 ierr = meshInstance->get_adjacencies ( &source[is], 1, 0, false, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 610 } 611 else if ( target_type == moab::MBEDGE ) 612 { 613 ierr = meshInstance->get_adjacencies ( &source[is], 1, 1, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 614 } 615 else if ( target_type <= moab::MBPOLYGON ) 616 { 617 ierr = meshInstance->get_adjacencies ( &source[is], 1, 2, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 618 } 619 else if ( target_type < moab::MBENTITYSET ) 620 { 621 ierr = meshInstance->get_adjacencies ( &source[is], 1, 3, true, adjents, moab::Interface::INTERSECT ); //MB_CHK_ERR_RET(ierr); 622 } 623 else /* Either EntitySet or MBMaxType -- Failures */ 624 { 625 MSQ_SETERR ( err ) ( process_itaps_error ( 1 ), MsqError::NOT_IMPLEMENTED ); // "Invalid Target entity type specified" 626 return; 627 } 628 629 ptr2[is + 1] = iadjoff + adjents.size(); 630 631 for ( unsigned iadj = 0; iadj < adjents.size(); ++iadj, ++iadjoff ) 632 { mArray.push_back ( adjents[iadj] ); } 633 } 634 635 num_adj = mArray.size(); 636 target.resize ( num_adj ); 637 std::copy ( mArray.begin(), mArray.end(), reinterpret_cast<moab::EntityHandle*> ( arrptr ( target ) ) ); 638 mArray.clear(); 639 } 640 641 if ( expand ) 642 { 643 for ( size_t i = num_offset; i > 0; --i ) 644 { offsets[i - 1] = static_cast<size_t> ( ptr2[i - 1] ); } 645 646 free ( ptr2 ); 647 } 648 649 // assert( (unsigned)num_offset == offsets.size() ); 650 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 651 } 652 653 vertices_get_attached_elements(const VertexHandle * vertices,size_t num_vertex,std::vector<ElementHandle> & elements,std::vector<size_t> & offsets,MsqError & err)654 void MsqMOAB::vertices_get_attached_elements ( 655 const VertexHandle* vertices, 656 size_t num_vertex, 657 std::vector<ElementHandle>& elements, 658 std::vector<size_t>& offsets, 659 MsqError& err ) 660 { 661 // moab::ErrorCode ierr; 662 bool cont; 663 assert ( sizeof ( EntityHandle ) == sizeof ( moab::EntityHandle ) ); 664 const moab::EntityHandle* verts = reinterpret_cast<const moab::EntityHandle*> ( vertices ); 665 get_adjacent_entities ( verts, num_vertex, inputSetType, elements, offsets, err ); 666 MSQ_ERRRTN ( err ); 667 668 moab::EntityHandle root_set = 0 /*meshInstance->get_root_set()*/; 669 670 // Remove all elements not in inputSet 671 if ( root_set != inputSet ) 672 { 673 std::vector<size_t>::iterator offset_iter = offsets.begin(); 674 size_t read_idx, write_idx; 675 moab::EntityHandle* bh = reinterpret_cast<moab::EntityHandle*> ( arrptr ( elements ) ); 676 677 for ( read_idx = write_idx = 0; read_idx < elements.size(); ++read_idx ) 678 { 679 if ( *offset_iter == read_idx ) 680 { 681 *offset_iter = write_idx; 682 ++offset_iter; 683 } 684 685 cont = meshInstance->contains_entities ( inputSet, &bh[read_idx], 1 ); 686 687 if ( cont ) 688 { elements[write_idx++] = elements[read_idx]; } 689 } 690 691 *offset_iter = write_idx; 692 elements.resize ( write_idx ); 693 } 694 } 695 696 697 //**************** Element Topology ***************** 698 699 700 /** Get connectivity 701 *\param elements - Array of length num_elems containing elements 702 * handles of elements for which connectivity is to 703 * be queried. 704 *\param vertices - Array of vertex handles in connectivity list. 705 *\param offsets - Indices into \ref vertex_handles, one per element 706 */ elements_get_attached_vertices(const ElementHandle * elements,size_t num_elems,std::vector<VertexHandle> & vertices,std::vector<size_t> & offsets,MBMesquite::MsqError & err)707 void MsqMOAB::elements_get_attached_vertices ( 708 const ElementHandle* elements, 709 size_t num_elems, 710 std::vector<VertexHandle>& vertices, 711 std::vector<size_t>& offsets, 712 MBMesquite::MsqError& err ) 713 { 714 moab::ErrorCode ierr; 715 assert ( sizeof ( moab::EntityHandle ) == sizeof ( EntityHandle ) ); 716 const moab::EntityHandle* elems = reinterpret_cast<const moab::EntityHandle*> ( elements ); 717 // get_adjacent_entities( elems, num_elems, moab::MBVERTEX, vertices, offsets, err ); MSQ_CHKERR(err); 718 offsets.resize ( num_elems + 1 ); 719 offsets[0] = 0; 720 vertices.clear(); 721 std::vector<moab::EntityHandle> mbverts; 722 723 for ( unsigned ie = 0; ie < num_elems; ++ie ) 724 { 725 const moab::EntityHandle* conn; 726 int nconn; 727 ierr = meshInstance->get_connectivity ( elems[ie], conn, nconn, false ); MB_CHK_ERR_RET ( ierr ); 728 offsets[ie + 1] = offsets[ie] + nconn; 729 730 for ( int iconn = 0; iconn < nconn; ++iconn ) 731 { 732 mbverts.push_back ( conn[iconn] ); 733 } 734 } 735 736 vertices.resize ( mbverts.size() ); 737 moab::EntityHandle* verts = reinterpret_cast<moab::EntityHandle*> ( arrptr ( vertices ) ); 738 739 for ( size_t iverts = 0; iverts < mbverts.size(); ++iverts ) 740 { verts[iverts] = mbverts[iverts]; } 741 742 mbverts.clear(); 743 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 744 } 745 746 get_all_elements(std::vector<ElementHandle> & elements,MsqError & err)747 void MsqMOAB::get_all_elements ( std::vector<ElementHandle>& elements, 748 MsqError& err ) 749 { 750 moab::ErrorCode ierr; 751 752 if ( inputSetType == moab::MBMAXTYPE ) 753 { 754 int num_vol, num_face; 755 756 ierr = meshInstance->get_number_entities_by_dimension ( inputSet, 2, num_face, false ); MB_CHK_ERR_RET ( ierr ); 757 ierr = meshInstance->get_number_entities_by_dimension ( inputSet, 3, num_vol, false ); MB_CHK_ERR_RET ( ierr ); 758 elements.resize ( num_face + num_vol ); 759 760 if ( elements.empty() ) 761 { return; } 762 763 moab::EntityHandle* ptr = reinterpret_cast<moab::EntityHandle*> ( arrptr ( elements ) ); 764 765 if ( num_face ) 766 { 767 std::vector<moab::EntityHandle> faces; 768 ierr = meshInstance->get_entities_by_dimension ( inputSet, 2, faces, false ); MB_CHK_ERR_RET ( ierr ); 769 assert ( faces.size() - num_face == 0 ); 770 std::copy ( faces.begin(), faces.end(), ptr ); 771 } 772 773 if ( num_vol ) 774 { 775 ptr += num_face; 776 std::vector<moab::EntityHandle> regions; 777 ierr = meshInstance->get_entities_by_dimension ( inputSet, 3, regions, false ); MB_CHK_ERR_RET ( ierr ); 778 assert ( regions.size() - num_vol == 0 ); 779 std::copy ( regions.begin(), regions.end(), ptr ); 780 } 781 } 782 else 783 { 784 int count; 785 ierr = meshInstance->get_number_entities_by_type ( inputSet, inputSetType, count, false ); MB_CHK_ERR_RET ( ierr ); 786 787 if ( !count ) 788 { return; } 789 790 elements.resize ( count ); 791 792 moab::EntityHandle* ptr = reinterpret_cast<moab::EntityHandle*> ( arrptr ( elements ) ); 793 std::vector<moab::EntityHandle> entities; 794 ierr = meshInstance->get_entities_by_type ( inputSet, inputSetType, entities, false ); MB_CHK_ERR_RET ( ierr ); 795 assert ( entities.size() - count == 0 ); 796 std::copy ( entities.begin(), entities.end(), ptr ); 797 } 798 799 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 800 } 801 get_all_vertices(std::vector<VertexHandle> & vertices,MsqError & err)802 void MsqMOAB::get_all_vertices ( std::vector<VertexHandle>& vertices, 803 MsqError& err ) 804 { 805 std::vector<ElementHandle> elems; 806 get_all_elements ( elems, err ); MSQ_CHKERR ( err ); 807 808 if ( elems.empty() ) 809 { return; } 810 811 std::vector<size_t> offsets; 812 elements_get_attached_vertices ( arrptr ( elems ), elems.size(), vertices, offsets, err ); MSQ_CHKERR ( err ); 813 814 std::sort ( vertices.begin(), vertices.end() ); 815 vertices.erase ( std::unique ( vertices.begin(), vertices.end() ), vertices.end() ); 816 } 817 818 819 // Returns the topologies of the given entities. The "entity_topologies" 820 // array must be at least "num_elements" in size. elements_get_topologies(const ElementHandle * element_handle_array,EntityTopology * element_topologies,size_t num_elements,MsqError & err)821 void MsqMOAB::elements_get_topologies ( 822 const ElementHandle* element_handle_array, 823 EntityTopology* element_topologies, 824 size_t num_elements, MsqError& err ) 825 { 826 if ( !num_elements ) 827 { return; } 828 829 assert ( sizeof ( ElementHandle ) == sizeof ( moab::EntityHandle ) ); 830 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( element_handle_array ); 831 832 for ( size_t i = 0; i < num_elements; ++i ) 833 { 834 element_topologies[i] = topologyMap[ meshInstance->type_from_handle ( arr[i] ) ]; 835 } 836 837 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 838 } 839 840 //**************** Memory Management **************** 841 // Tells the mesh that the client is finished with a given 842 // entity handle. release_entity_handles(const MBMesquite::Mesh::EntityHandle *,size_t,MsqError &)843 void MsqMOAB::release_entity_handles ( 844 const MBMesquite::Mesh::EntityHandle* /*handle_array*/, 845 size_t /*num_handles*/, MsqError& /*err*/ ) 846 { 847 // Do nothing 848 } 849 850 // Instead of deleting a Mesh when you think you are done, 851 // call release(). In simple cases, the implementation could 852 // just call the destructor. More sophisticated implementations 853 // may want to keep the Mesh object to live longer than Mesquite 854 // is using it. release()855 void MsqMOAB::release() 856 { 857 } 858 859 //**************** Tags **************** tag_create(const std::string & name,TagType type,unsigned length,const void *,MsqError & err)860 TagHandle MsqMOAB::tag_create ( const std::string& name, 861 TagType type, unsigned length, 862 const void*, 863 MsqError& err ) 864 { 865 moab::DataType moab_type; 866 867 switch ( type ) 868 { 869 case MBMesquite::Mesh::BYTE: moab_type = moab::MB_TYPE_OPAQUE; break; 870 871 case MBMesquite::Mesh::INT: moab_type = moab::MB_TYPE_INTEGER; break; 872 873 case MBMesquite::Mesh::DOUBLE: moab_type = moab::MB_TYPE_DOUBLE; break; 874 875 case MBMesquite::Mesh::HANDLE: moab_type = moab::MB_TYPE_HANDLE; break; 876 877 default: 878 MSQ_SETERR ( err ) ( "Invalid tag type", MsqError::INVALID_ARG ); 879 return 0; 880 } 881 882 moab::ErrorCode ierr; 883 moab::Tag result = 0; 884 ierr = meshInstance->tag_get_handle ( name.c_str(), length, moab_type, result, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE ); MB_CHK_ERR_RET_VAL ( ierr, result ); 885 886 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 887 return static_cast<TagHandle> ( result ); 888 } 889 tag_delete(TagHandle handle,MsqError & err)890 void MsqMOAB::tag_delete ( TagHandle handle, MsqError& err ) 891 { 892 moab::ErrorCode ierr = meshInstance->tag_delete ( static_cast<moab::Tag> ( handle ) ); MB_CHK_ERR_RET ( ierr ); 893 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 894 } 895 tag_get(const std::string & name,MsqError & err)896 TagHandle MsqMOAB::tag_get ( const std::string& name, MsqError& err ) 897 { 898 moab::ErrorCode ierr; 899 moab::Tag handle = 0; 900 ierr = meshInstance->tag_get_handle ( name.c_str(), handle ); MB_CHK_ERR_RET_VAL ( ierr, handle ); 901 902 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 903 return static_cast<TagHandle> ( handle ); 904 } 905 906 tag_properties(TagHandle handle,std::string & name_out,TagType & type_out,unsigned & length_out,MsqError & err)907 void MsqMOAB::tag_properties ( TagHandle handle, 908 std::string& name_out, 909 TagType& type_out, 910 unsigned& length_out, 911 MsqError& err ) 912 { 913 std::string buffer; 914 moab::ErrorCode ierr; 915 moab::DataType itype; 916 int size; 917 918 moab::Tag th = static_cast<moab::Tag> ( handle ); 919 ierr = meshInstance->tag_get_name ( th, buffer ); MB_CHK_ERR_RET ( ierr ); 920 ierr = meshInstance->tag_get_length ( th, size ); MB_CHK_ERR_RET ( ierr ); 921 ierr = meshInstance->tag_get_data_type ( th, itype ); MB_CHK_ERR_RET ( ierr ); 922 923 name_out = buffer; 924 length_out = size; 925 926 switch ( itype ) 927 { 928 case moab::MB_TYPE_OPAQUE : type_out = MBMesquite::Mesh::BYTE ; break; 929 930 case moab::MB_TYPE_INTEGER : type_out = MBMesquite::Mesh::INT ; break; 931 932 case moab::MB_TYPE_DOUBLE : type_out = MBMesquite::Mesh::DOUBLE; break; 933 934 case moab::MB_TYPE_HANDLE : type_out = MBMesquite::Mesh::HANDLE; break; 935 936 default: 937 MSQ_SETERR ( err ) ( "Unsupported iMesh tag type", MsqError::NOT_IMPLEMENTED ); 938 return; 939 } 940 941 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 942 } 943 tag_set_element_data(TagHandle tag,size_t num_elems,const ElementHandle * array,const void * data,MsqError & err)944 void MsqMOAB::tag_set_element_data ( TagHandle tag, 945 size_t num_elems, 946 const ElementHandle* array, 947 const void* data, 948 MsqError& err ) 949 { 950 tag_set_data ( tag, num_elems, array, data, err ); 951 } 952 tag_set_vertex_data(TagHandle tag,size_t num_elems,const VertexHandle * array,const void * data,MsqError & err)953 void MsqMOAB::tag_set_vertex_data ( TagHandle tag, 954 size_t num_elems, 955 const VertexHandle* array, 956 const void* data, 957 MsqError& err ) 958 { 959 tag_set_data ( tag, num_elems, array, data, err ); 960 } 961 tag_set_data(TagHandle tag,size_t num_elems,const EntityHandle * array,const void * data,MsqError & err)962 void MsqMOAB::tag_set_data ( TagHandle tag, 963 size_t num_elems, 964 const EntityHandle* array, 965 const void* data, 966 MsqError& err ) 967 { 968 moab::ErrorCode ierr; 969 int size; 970 moab::Tag th = static_cast<moab::Tag> ( tag ); 971 ierr = meshInstance->tag_get_length ( th, size ); MB_CHK_ERR_RET ( ierr ); 972 973 assert ( sizeof ( EntityHandle ) == sizeof ( moab::EntityHandle ) ); 974 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( array ); 975 ierr = meshInstance->tag_set_data ( th, arr, num_elems, data ); MB_CHK_ERR_RET ( ierr ); 976 977 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 978 } 979 980 tag_get_element_data(TagHandle tag,size_t num_elems,const ElementHandle * array,void * data,MsqError & err)981 void MsqMOAB::tag_get_element_data ( TagHandle tag, 982 size_t num_elems, 983 const ElementHandle* array, 984 void* data, 985 MsqError& err ) 986 { 987 tag_get_data ( tag, num_elems, array, data, err ); 988 } 989 tag_get_vertex_data(TagHandle tag,size_t num_verts,const VertexHandle * array,void * data,MsqError & err)990 void MsqMOAB::tag_get_vertex_data ( TagHandle tag, 991 size_t num_verts, 992 const VertexHandle* array, 993 void* data, 994 MsqError& err ) 995 { 996 tag_get_data ( tag, num_verts, array, data, err ); 997 } 998 tag_get_data(TagHandle tag,size_t num_elems,const EntityHandle * array,void * data,MsqError & err)999 void MsqMOAB::tag_get_data ( TagHandle tag, 1000 size_t num_elems, 1001 const EntityHandle* array, 1002 void* data, 1003 MsqError& err ) 1004 { 1005 moab::ErrorCode ierr; 1006 1007 #if IMESH_VERSION_ATLEAST(1,1) 1008 void* ptr = data; 1009 #else 1010 char* ptr = static_cast<char*> ( data ); 1011 #endif 1012 1013 assert ( sizeof ( EntityHandle ) == sizeof ( moab::EntityHandle ) ); 1014 moab::Tag th = static_cast<moab::Tag> ( tag ); 1015 const moab::EntityHandle* arr = reinterpret_cast<const moab::EntityHandle*> ( array ); 1016 ierr = meshInstance->tag_get_data ( th, arr, num_elems, ptr ); MB_CHK_ERR_RET ( ierr ); 1017 1018 err.set_error ( MBMesquite::MsqError::NO_ERROR ); 1019 } 1020 1021 } // namespace Mesquite 1022 1023