moab
|
00001 #include "iMeshP.h" 00002 #include "moab_mpi.h" 00003 #include <iostream> 00004 #include <algorithm> 00005 #include <vector> 00006 #include <sstream> 00007 #include <assert.h> 00008 #include <math.h> 00009 #include <map> 00010 #include <string.h> 00011 #include <stdio.h> // remove() 00012 00013 #if !defined(_MSC_VER) && !defined(__MINGW32__) 00014 #include <unistd.h> 00015 #endif 00016 00017 #define STRINGIFY_(X) #X 00018 #define STRINGIFY(X) STRINGIFY_(X) 00019 const char* const FILENAME = "iMeshP_test_file"; 00020 00021 00022 /************************************************************************** 00023 Error Checking 00024 **************************************************************************/ 00025 00026 #define CHKERR do { \ 00027 if (ierr) { \ 00028 std::cerr << "Error code " << ierr << " at " << __FILE__ << ":" << __LINE__ << std::endl;\ 00029 return ierr; \ 00030 } \ 00031 } while (false) 00032 00033 #define PCHECK do { ierr = is_any_proc_error(ierr); CHKERR; } while(false) 00034 00035 // Use my_rank instead of rank to avoid shadowed declaration 00036 #define ASSERT(A) do { \ 00037 if (is_any_proc_error(!(A))) { \ 00038 int my_rank = 0; \ 00039 MPI_Comm_rank( MPI_COMM_WORLD, &my_rank ); \ 00040 if (0 == my_rank) std::cerr << "Failed assertion: " #A << std::endl \ 00041 << " at " __FILE__ ":" << __LINE__ << std::endl; \ 00042 return 1; \ 00043 } } while (false) 00044 00045 // Test if is_my_error is non-zero on any processor in MPI_COMM_WORLD 00046 int is_any_proc_error( int is_my_error ) 00047 { 00048 int result; 00049 int err = MPI_Allreduce( &is_my_error, &result, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD ); 00050 return err || result; 00051 } 00052 00053 /************************************************************************** 00054 Test Declarations 00055 **************************************************************************/ 00056 00057 class PartMap; 00058 00063 int test_load( iMesh_Instance, iMeshP_PartitionHandle prtn, PartMap& map, int comm_size ); 00064 00065 00073 int test_get_partitions( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00074 00082 int test_get_parts( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00083 00092 int test_get_by_type( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00093 00102 int test_get_by_topo( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00103 00114 int test_part_id_handle( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00115 00122 int test_part_rank( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00123 00132 int test_get_neighbors( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00133 00140 int test_get_part_boundary( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00141 00148 int test_part_boundary_iter( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00149 00155 int test_get_adjacencies( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00156 00163 int test_entity_iterator( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00164 00173 int test_entity_owner( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00174 00181 int test_entity_status( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00182 00189 int test_entity_copy_parts( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00190 00198 int test_entity_copies( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00199 00205 int test_create_ghost_ents( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00206 00213 int test_push_tag_data_iface( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00214 int test_push_tag_data_ghost( iMesh_Instance, iMeshP_PartitionHandle prtn, const PartMap& ); 00215 00216 int test_exchange_ents( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ); 00217 00218 /************************************************************************** 00219 Helper Funcions 00220 **************************************************************************/ 00221 00222 class PartMap 00223 { 00224 public: 00225 int num_parts() const 00226 { return sortedPartList.size(); } 00227 00228 iMeshP_Part part_id_from_local_id( int local_id ) const 00229 { return sortedPartList[idx_from_local_id(local_id)]; } 00230 00231 int local_id_from_part_id( iMeshP_Part part ) const 00232 { return partLocalIds[idx_from_part_id(part)]; } 00233 00234 int rank_from_part_id( iMeshP_Part part ) const 00235 { return partRanks[idx_from_part_id(part)]; } 00236 00237 int rank_from_local_id( int id ) const 00238 { return partRanks[idx_from_local_id(id)]; } 00239 00240 int count_from_rank( int rank ) const 00241 { return std::count( partRanks.begin(), partRanks.end(), rank ); } 00242 00243 void part_id_from_rank( int rank, std::vector<iMeshP_Part>& parts ) const; 00244 00245 void local_id_from_rank( int rank, std::vector<int>& ids ) const; 00246 00247 const std::vector<iMeshP_Part>& get_parts() const 00248 { return sortedPartList; } 00249 00250 const std::vector<int>& get_ranks() const 00251 { return partRanks; } 00252 00253 int build_map( iMesh_Instance imesh, 00254 iMeshP_PartitionHandle partition, 00255 int num_expected_parts ); 00256 00257 static int part_from_coords( iMesh_Instance imesh, 00258 iMeshP_PartHandle part, 00259 int& id_out ); 00260 00261 private: 00262 inline int idx_from_part_id( iMeshP_Part id ) const 00263 { return std::lower_bound( sortedPartList.begin(), sortedPartList.end(), id ) 00264 - sortedPartList.begin(); } 00265 inline int idx_from_local_id( int id ) const 00266 { return localIdReverseMap[id]; } 00267 00268 std::vector<iMeshP_Part> sortedPartList; 00269 std::vector<int> partRanks; 00270 std::vector<int> partLocalIds; 00271 std::vector<int> localIdReverseMap; 00272 }; 00273 00275 int create_mesh( const char* filename, int num_parts ); 00276 00278 int vertex_tag( iMesh_Instance imesh, iBase_EntityHandle vertex, int& tag ); 00279 00280 int get_local_parts( iMesh_Instance instance, 00281 iMeshP_PartitionHandle prtn, 00282 std::vector<iMeshP_PartHandle>& handles, 00283 std::vector<iMeshP_Part>* ids = 0 ) 00284 { 00285 iMeshP_PartHandle* arr = 0; 00286 int ierr, alloc = 0, size; 00287 iMeshP_getLocalParts( instance, prtn, &arr, &alloc, &size, &ierr ); 00288 CHKERR; 00289 handles.resize( size ); 00290 std::copy( arr, arr + size, handles.begin() ); 00291 if (!ids) 00292 return iBase_SUCCESS; 00293 00294 ids->resize( size ); 00295 alloc = size; 00296 iMeshP_Part* ptr = &(*ids)[0]; 00297 iMeshP_getPartIdsFromPartHandlesArr( instance, prtn, &handles[0], handles.size(), 00298 &ptr, &alloc, &size, &ierr ); 00299 CHKERR; 00300 assert( size == (int)ids->size() ); 00301 assert( ptr == &(*ids)[0] ); 00302 return iBase_SUCCESS; 00303 } 00304 00305 00306 static int get_entities( iMesh_Instance imesh, 00307 iBase_EntitySetHandle set, 00308 iBase_EntityType type, 00309 iMesh_EntityTopology topo, 00310 std::vector<iBase_EntityHandle>& entities ) 00311 { 00312 iBase_EntityHandle* array = 0; 00313 int junk = 0, size = 0, err; 00314 iMesh_getEntities( imesh, set, type, topo, &array, &junk, &size, &err ); 00315 if (!err) { 00316 entities.clear(); 00317 entities.resize( size ); 00318 std::copy( array, array + size, entities.begin() ); 00319 free( array ); 00320 } 00321 return err; 00322 } 00323 00324 static int get_part_quads_and_verts( iMesh_Instance imesh, 00325 iMeshP_PartHandle part, 00326 std::vector<iBase_EntityHandle>& elems, 00327 std::vector<iBase_EntityHandle>& verts ) 00328 { 00329 int ierr = get_entities( imesh, part, iBase_FACE, iMesh_QUADRILATERAL, elems ); 00330 CHKERR; 00331 00332 verts.resize(4*elems.size()); 00333 std::vector<int> junk(elems.size()+1); 00334 int junk1 = verts.size(), count, junk2 = junk.size(), junk3; 00335 iBase_EntityHandle* junk4 = &verts[0]; 00336 int* junk5 = &junk[0]; 00337 iMesh_getEntArrAdj( imesh, &elems[0], elems.size(), iBase_VERTEX, 00338 &junk4, &junk1, &count, 00339 &junk5, &junk2, &junk3, &ierr ); 00340 CHKERR; 00341 assert( junk1 == (int)verts.size() ); 00342 assert( count == (int)(4*elems.size()) ); 00343 assert( junk2 == (int)junk.size() ); 00344 assert( junk4 == &verts[0] ); 00345 assert( junk5 == &junk[0] ); 00346 std::sort( verts.begin(), verts.end() ); 00347 verts.erase( std::unique( verts.begin(), verts.end() ), verts.end() ); 00348 return iBase_SUCCESS; 00349 } 00350 00351 static int get_coords( iMesh_Instance imesh, 00352 const iBase_EntityHandle* verts, 00353 int num_verts, 00354 double* coords ) 00355 { 00356 double* junk1 = coords; 00357 int junk2 = 3*num_verts; 00358 int junk3; 00359 int ierr; 00360 iMesh_getVtxArrCoords( imesh, verts, num_verts, iBase_INTERLEAVED, &junk1, &junk2, &junk3, &ierr ); 00361 if (iBase_SUCCESS != ierr) 00362 return ierr; 00363 assert( junk1 == coords ); 00364 assert( junk2 == 3*num_verts ); 00365 assert( junk3 == 3*num_verts ); 00366 return iBase_SUCCESS; 00367 } 00368 00369 /************************************************************************** 00370 Main Method 00371 **************************************************************************/ 00372 00373 #define RUN_TEST(A) run_test( &A, #A ) 00374 00375 int run_test( int (*func)(iMesh_Instance, iMeshP_PartitionHandle, const PartMap&), 00376 const char* func_name ) 00377 { 00378 int rank, size, ierr; 00379 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00380 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00381 00382 if (rank == 0) { 00383 ierr = create_mesh( FILENAME, size ); 00384 } 00385 MPI_Bcast( &ierr, 1, MPI_INT, 0, MPI_COMM_WORLD ); 00386 if (ierr) { 00387 if (rank == 0) { 00388 std::cerr << "Failed to create input test file on root processor. Aborting." 00389 << std::endl; 00390 } 00391 abort(); 00392 } 00393 00394 iMesh_Instance imesh; 00395 iMesh_newMesh( 0, &imesh, &ierr, 0 ); 00396 PCHECK; 00397 00398 iMeshP_PartitionHandle prtn; 00399 iMeshP_createPartitionAll( imesh, MPI_COMM_WORLD, &prtn, &ierr ); 00400 PCHECK; 00401 00402 PartMap map; 00403 ierr = test_load( imesh, prtn, map, size ); 00404 if (ierr) { 00405 if (rank == 0) { 00406 std::cerr << "Failed to load input mesh." << std::endl 00407 << "Cannot run further tests." << std::endl 00408 << "ABORTING" << std::endl; 00409 } 00410 abort(); 00411 } 00412 00413 int result = (*func)(imesh,prtn,map); 00414 int is_err = is_any_proc_error( result ); 00415 if (rank == 0) { 00416 if (is_err) 00417 std::cout << func_name << " : FAILED!!" << std::endl; 00418 else 00419 std::cout << func_name << " : success" << std::endl; 00420 } 00421 00422 iMesh_dtor( imesh, &ierr ); 00423 CHKERR; 00424 return is_err; 00425 } 00426 00427 int main( int argc, char* argv[] ) 00428 { 00429 MPI_Init(&argc, &argv); 00430 int size, rank; 00431 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00432 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00433 00434 if (argc > 2 && !strcmp(argv[1], "-p")) { 00435 #if !defined(_MSC_VER) && !defined(__MINGW32__) 00436 std::cout << "Processor " << rank << " of " << size << " with PID " << getpid() << std::endl; 00437 std::cout.flush(); 00438 #endif 00439 // loop forever on requested processor, giving the user time 00440 // to attach a debugger. Once the debugger in attached, user 00441 // can change 'pause'. E.g. on gdb do "set var pause = 0" 00442 if (atoi(argv[2]) == rank) { 00443 volatile int pause = 1; 00444 while (pause); 00445 } 00446 MPI_Barrier(MPI_COMM_WORLD); 00447 } 00448 00449 int num_errors = 0; 00450 num_errors += RUN_TEST( test_get_partitions ); 00451 num_errors += RUN_TEST( test_get_parts ); 00452 num_errors += RUN_TEST( test_get_by_type ); 00453 num_errors += RUN_TEST( test_get_by_topo ); 00454 num_errors += RUN_TEST( test_part_id_handle ); 00455 num_errors += RUN_TEST( test_part_rank ); 00456 num_errors += RUN_TEST( test_get_neighbors ); 00457 num_errors += RUN_TEST( test_get_part_boundary ); 00458 num_errors += RUN_TEST( test_part_boundary_iter ); 00459 num_errors += RUN_TEST( test_get_adjacencies ); 00460 num_errors += RUN_TEST( test_entity_iterator ); 00461 num_errors += RUN_TEST( test_entity_owner ); 00462 num_errors += RUN_TEST( test_entity_status ); 00463 num_errors += RUN_TEST( test_entity_copy_parts ); 00464 num_errors += RUN_TEST( test_entity_copies ); 00465 num_errors += RUN_TEST( test_push_tag_data_iface ); 00466 num_errors += RUN_TEST( test_push_tag_data_ghost ); 00467 num_errors += RUN_TEST( test_create_ghost_ents ); 00468 num_errors += RUN_TEST( test_exchange_ents ); 00469 00470 // wait until all procs are done before writing summary data 00471 std::cout.flush(); 00472 MPI_Barrier( MPI_COMM_WORLD ); 00473 00474 // clean up output file 00475 if (rank == 0) 00476 remove( FILENAME ); 00477 00478 if (rank == 0) { 00479 if (!num_errors) 00480 std::cout << "All tests passed" << std::endl; 00481 else 00482 std::cout << num_errors << " TESTS FAILED!" << std::endl; 00483 } 00484 00485 MPI_Finalize(); 00486 return num_errors; 00487 } 00488 00489 // Create a mesh 00490 // 00491 // 00492 // Groups of four quads will be arranged into parts as follows: 00493 // +------+------+------+------+------+----- 00494 // | | | 00495 // | | | 00496 // + Part 0 + Part 2 + Part 4 00497 // | | | 00498 // | | | 00499 // +------+------+------+------+------+----- 00500 // | | | 00501 // | | | 00502 // + Part 1 + Part 3 + Part 5 00503 // | | | 00504 // | | | 00505 // +------+------+------+------+------+----- 00506 // 00507 // Vertices will be enumerated as follows: 00508 // 1------6-----11-----16-----21-----26----- 00509 // | | | 00510 // | | | 00511 // 2 7 12 17 22 27 00512 // | | | 00513 // | | | 00514 // 3------8-----13-----18-----23-----28----- 00515 // | | | 00516 // | | | 00517 // 4 9 14 19 24 29 00518 // | | | 00519 // | | | 00520 // 5-----10-----15-----20-----25-----30----- 00521 // 00522 // Element IDs will be [4*rank+1,4*rank+5] 00523 template <int size> struct EHARR 00524 { 00525 iBase_EntityHandle h[size]; 00526 iBase_EntityHandle& operator[](int i){ return h[i]; } 00527 operator iBase_EntityHandle*() { return h; } 00528 }; 00529 int create_mesh( const char* filename, int num_parts ) 00530 { 00531 const char* tagname = "GLOBAL_ID"; 00532 int ierr; 00533 00534 iMesh_Instance imesh; 00535 iMesh_newMesh( 0, &imesh, &ierr, 0 ); CHKERR; 00536 00537 const int num_full_cols = 2 * (num_parts / 2); 00538 const int need_half_cols = num_parts % 2; 00539 const int num_cols = num_full_cols + 2*need_half_cols; 00540 const int num_vtx = 5+5*num_cols - 4*(num_parts%2); 00541 std::vector< EHARR<5> > vertices( num_cols + 1 ); 00542 std::vector< EHARR<4> > elements( num_cols ); 00543 std::vector<int> vertex_ids( num_vtx ); 00544 std::vector<iBase_EntityHandle> vertex_list(num_vtx); 00545 for (int i = 0; i < num_vtx; ++i) 00546 vertex_ids[i] = i+1; 00547 00548 // create vertices 00549 int vl_pos = 0; 00550 for (int i = 0; i <= num_cols; ++i) { 00551 double coords[15] = { static_cast<double>(i), 0, 0, 00552 static_cast<double>(i), 1, 0, 00553 static_cast<double>(i), 2, 0, 00554 static_cast<double>(i), 3, 0, 00555 static_cast<double>(i), 4, 0 }; 00556 iBase_EntityHandle* ptr = vertices[i]; 00557 const int n = (num_full_cols == num_cols || i <= num_full_cols) ? 5 : 3; 00558 int junk1 = n, junk2 = n; 00559 iMesh_createVtxArr( imesh, n, iBase_INTERLEAVED, coords, 3*n, 00560 &ptr, &junk1, &junk2, &ierr ); CHKERR; 00561 assert( ptr == vertices[i] ); 00562 assert( junk1 == n ); 00563 assert( junk2 == n ); 00564 for (int j = 0; j < n; ++j) 00565 vertex_list[vl_pos++] = vertices[i][j]; 00566 } 00567 00568 // create elements 00569 for (int i = 0; i < num_cols; ++i) { 00570 iBase_EntityHandle conn[16]; 00571 for (int j = 0; j < 4; ++j) { 00572 conn[4*j ] = vertices[i ][j ]; 00573 conn[4*j+1] = vertices[i ][j+1]; 00574 conn[4*j+2] = vertices[i+1][j+1]; 00575 conn[4*j+3] = vertices[i+1][j ]; 00576 } 00577 iBase_EntityHandle* ptr = elements[i]; 00578 const int n = (i < num_full_cols) ? 4 : 2; 00579 int junk1 = n, junk2 = n, junk3 = n, junk4 = n; 00580 int stat[4]; 00581 int* ptr2 = stat; 00582 iMesh_createEntArr( imesh, 00583 iMesh_QUADRILATERAL, 00584 conn, 4*n, 00585 &ptr, &junk1, &junk2, 00586 &ptr2, &junk3, &junk4, 00587 &ierr ); CHKERR; 00588 assert( ptr == elements[i] ); 00589 assert( junk1 == n ); 00590 assert( junk2 == n ); 00591 assert( ptr2 == stat ); 00592 assert( junk3 == n ); 00593 assert( junk4 == n ); 00594 } 00595 00596 // create partition 00597 iMeshP_PartitionHandle partition; 00598 iMeshP_createPartitionAll( imesh, MPI_COMM_SELF, &partition, &ierr ); CHKERR; 00599 for (int i = 0; i < num_parts; ++i) { 00600 iMeshP_PartHandle part; 00601 iMeshP_createPart( imesh, partition, &part, &ierr ); CHKERR; 00602 iBase_EntityHandle quads[] = { elements[2*(i/2) ][2*(i%2) ], 00603 elements[2*(i/2)+1][2*(i%2) ], 00604 elements[2*(i/2) ][2*(i%2)+1], 00605 elements[2*(i/2)+1][2*(i%2)+1] }; 00606 iMesh_addEntArrToSet( imesh, quads, 4, part, &ierr ); CHKERR; 00607 } 00608 00609 // assign global ids to vertices 00610 iBase_TagHandle id_tag = 0; 00611 iMesh_getTagHandle( imesh, tagname, &id_tag, &ierr, strlen(tagname) ); 00612 if (iBase_SUCCESS == ierr) { 00613 int tag_size, tag_type; 00614 iMesh_getTagSizeValues( imesh, id_tag, &tag_size, &ierr ); 00615 CHKERR; 00616 if (tag_size != 1) 00617 return iBase_TAG_ALREADY_EXISTS; 00618 iMesh_getTagType( imesh, id_tag, &tag_type, &ierr ); 00619 CHKERR; 00620 if (tag_type != iBase_INTEGER) 00621 return iBase_TAG_ALREADY_EXISTS; 00622 } 00623 else { 00624 iMesh_createTag( imesh, tagname, 1, iBase_INTEGER, &id_tag, &ierr, strlen(tagname) ); 00625 CHKERR; 00626 } 00627 iMesh_setIntArrData( imesh, &vertex_list[0], num_vtx, id_tag, &vertex_ids[0], num_vtx, &ierr ); 00628 CHKERR; 00629 00630 // write file 00631 iBase_EntitySetHandle root_set; 00632 iMesh_getRootSet( imesh, &root_set, &ierr ); 00633 iMeshP_saveAll( imesh, partition, root_set, filename, 0, &ierr, strlen(filename), 0 ); 00634 CHKERR; 00635 00636 iMesh_dtor( imesh, &ierr ); CHKERR; 00637 00638 return 0; 00639 } 00640 00641 // generate unique for each vertex from coordinates. 00642 // Assume integer coordinate values with x in [0,inf] and y in [0,4] 00643 // as generated by create_mean(..). 00644 int vertex_tag( iMesh_Instance imesh, iBase_EntityHandle vertex, int& tag ) 00645 { 00646 int ierr; 00647 double x, y, z; 00648 iMesh_getVtxCoord( imesh, vertex, &x, &y, &z, &ierr ); 00649 CHKERR; 00650 00651 int xc = (int)round(x); 00652 int yc = (int)round(y); 00653 tag = 5*xc + yc + 1; 00654 return ierr; 00655 } 00656 00657 /************************************************************************** 00658 Test Implementations 00659 **************************************************************************/ 00660 00661 int test_load( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, PartMap& map, int proc_size ) 00662 { 00663 int ierr; 00664 00665 iBase_EntitySetHandle root_set; 00666 iMesh_getRootSet( imesh, &root_set, &ierr ); 00667 const char *opt = "moab:PARTITION=PARALLEL_PARTITION"; 00668 iMeshP_loadAll( imesh, prtn, root_set, FILENAME, opt, &ierr, strlen(FILENAME), strlen(opt) ); 00669 PCHECK; 00670 00671 00672 ierr = map.build_map( imesh, prtn, proc_size ); 00673 CHKERR; 00674 return iBase_SUCCESS; 00675 } 00676 00677 00685 int test_get_partitions( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& ) 00686 { 00687 int ierr; 00688 00689 // test iMeshP_getPartitionCom 00690 MPI_Comm comm = MPI_COMM_SELF; 00691 iMeshP_getPartitionComm( imesh, prtn, &comm, &ierr ); 00692 PCHECK; 00693 ASSERT(comm == MPI_COMM_WORLD); 00694 00695 00696 // test iMeshP_getPartitions 00697 iMeshP_PartitionHandle* array = 0; 00698 int alloc = 0, size = -1; 00699 iMeshP_getPartitions( imesh, &array, &alloc, &size, &ierr ); 00700 PCHECK; 00701 ASSERT(array != 0); 00702 ASSERT(alloc == size); 00703 ASSERT(size > 0); 00704 int idx = std::find(array, array+size, prtn) - array; 00705 free(array); 00706 ASSERT(idx < size); 00707 00708 // test iMesP_getNumPartitions 00709 int size2 = -1; 00710 iMeshP_getNumPartitions( imesh, &size2, &ierr ); 00711 PCHECK; 00712 ASSERT(size2 == size); 00713 return 0; 00714 } 00715 00716 00717 00725 int test_get_parts( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 00726 { 00727 int size, rank, ierr; 00728 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00729 MPI_Comm_size( MPI_COMM_WORLD, &size ); 00730 00731 int num_part_g; 00732 iMeshP_getNumGlobalParts( imesh, prtn, &num_part_g, &ierr ); 00733 PCHECK; 00734 ASSERT( num_part_g == map.num_parts() ); 00735 00736 int num_part_l; 00737 iMeshP_getNumLocalParts( imesh, prtn, &num_part_l, &ierr ); 00738 PCHECK; 00739 ASSERT( num_part_l == map.count_from_rank( rank ) ); 00740 00741 std::vector<iMeshP_PartHandle> parts(num_part_l); 00742 iMeshP_PartHandle* ptr = &parts[0]; 00743 int junk1 = num_part_l, count = -1; 00744 iMeshP_getLocalParts( imesh, prtn, &ptr, &junk1, &count, &ierr ); 00745 PCHECK; 00746 assert( ptr == &parts[0] ); 00747 assert( junk1 == num_part_l ); 00748 ASSERT( count == num_part_l ); 00749 00750 return iBase_SUCCESS; 00751 } 00752 00753 static int test_get_by_type_topo_all( iMesh_Instance imesh, 00754 iMeshP_PartitionHandle prtn, 00755 bool test_type, 00756 int num_parts ) 00757 { 00758 // calculate number of quads and vertices in entire mesh 00759 // from number of parts (see create_mesh(..) function.) 00760 const int expected_global_quad_count = 4 * num_parts; 00761 const int num_col = 2 * (num_parts / 2 + num_parts % 2); 00762 const int expected_global_vtx_count = num_parts == 1 ? 9 : 00763 num_parts % 2 ? 1 + 5*num_col : 00764 5 + 5*num_col; 00765 00766 // test getNumOf*All for root set 00767 int ierr, count; 00768 iBase_EntitySetHandle root; 00769 iMesh_getRootSet( imesh, &root, &ierr ); 00770 if (test_type) 00771 iMeshP_getNumOfTypeAll( imesh, prtn, root, iBase_VERTEX, &count, &ierr ); 00772 else 00773 iMeshP_getNumOfTopoAll( imesh, prtn, root, iMesh_POINT, &count, &ierr ); 00774 PCHECK; 00775 ASSERT( count == expected_global_vtx_count ); 00776 if (test_type) 00777 iMeshP_getNumOfTypeAll( imesh, prtn, root, iBase_FACE, &count, &ierr ); 00778 else 00779 iMeshP_getNumOfTopoAll( imesh, prtn, root, iMesh_QUADRILATERAL, &count, &ierr ); 00780 PCHECK; 00781 ASSERT( count == expected_global_quad_count ); 00782 00783 // create an entity set containing half of the quads 00784 std::vector<iBase_EntityHandle> all_quads, half_quads; 00785 ierr = get_entities( imesh, root, iBase_FACE, iMesh_QUADRILATERAL, all_quads ); 00786 assert( 0 == all_quads.size() % 2 ); 00787 half_quads.resize(all_quads.size()/2); 00788 for (size_t i = 0; i < all_quads.size() / 2; ++i) 00789 half_quads[i] = all_quads[2*i]; 00790 iBase_EntitySetHandle set; 00791 iMesh_createEntSet( imesh, 1, &set, &ierr ); 00792 CHKERR; 00793 iMesh_addEntArrToSet( imesh, &half_quads[0], half_quads.size(), set, &ierr ); 00794 CHKERR; 00795 00796 // test getNumOf*All with defined set 00797 if (test_type) 00798 iMeshP_getNumOfTypeAll( imesh, prtn, set, iBase_VERTEX, &count, &ierr ); 00799 else 00800 iMeshP_getNumOfTopoAll( imesh, prtn, set, iMesh_POINT, &count, &ierr ); 00801 PCHECK; 00802 ASSERT( count == 0 ); 00803 if (test_type) 00804 iMeshP_getNumOfTypeAll( imesh, prtn, set, iBase_FACE, &count, &ierr ); 00805 else 00806 iMeshP_getNumOfTopoAll( imesh, prtn, set, iMesh_QUADRILATERAL, &count, &ierr ); 00807 PCHECK; 00808 ASSERT( count == expected_global_quad_count/2 ); 00809 00810 return 0; 00811 } 00812 00813 static int test_get_by_type_topo_local( iMesh_Instance imesh, 00814 iMeshP_PartitionHandle prtn, 00815 bool test_type ) 00816 { 00817 int ierr; 00818 iBase_EntitySetHandle root; 00819 iMesh_getRootSet( imesh, &root, &ierr ); 00820 00821 // select a single part 00822 std::vector<iMeshP_PartHandle> parts; 00823 ierr = get_local_parts( imesh, prtn, parts ); 00824 CHKERR; 00825 iMeshP_PartHandle part = parts.front(); 00826 00827 // get the entities contained in the part 00828 std::vector<iBase_EntityHandle> part_quads, part_all; 00829 ierr = get_entities( imesh, part, iBase_FACE, iMesh_QUADRILATERAL, part_quads ); CHKERR; 00830 ierr = get_entities( imesh, part, iBase_ALL_TYPES, iMesh_ALL_TOPOLOGIES, part_all ); CHKERR; 00831 00832 // compare local counts (using root set) 00833 00834 int count; 00835 if (test_type) 00836 iMeshP_getNumOfType( imesh, prtn, part, root, iBase_FACE, &count, &ierr ); 00837 else 00838 iMeshP_getNumOfTopo( imesh, prtn, part, root, iMesh_QUADRILATERAL, &count, &ierr ); 00839 CHKERR; 00840 ASSERT( count == (int)part_quads.size() ); 00841 00842 if (test_type) 00843 iMeshP_getNumOfType( imesh, prtn, part, root, iBase_ALL_TYPES, &count, &ierr ); 00844 else 00845 iMeshP_getNumOfTopo( imesh, prtn, part, root, iMesh_ALL_TOPOLOGIES, &count, &ierr ); 00846 CHKERR; 00847 ASSERT( count == (int)part_all.size() ); 00848 00849 // compare local contents (using root set) 00850 00851 iBase_EntityHandle* ptr = 0; 00852 int num_ent, junk1 = 0; 00853 iMeshP_getEntities( imesh, prtn, part, root, test_type ? iBase_FACE : iBase_ALL_TYPES, 00854 test_type ? iMesh_ALL_TOPOLOGIES : iMesh_QUADRILATERAL, 00855 &ptr, &junk1, &num_ent, &ierr ); CHKERR; 00856 std::vector<iBase_EntityHandle> act_quads( ptr, ptr+num_ent ); 00857 free(ptr); 00858 junk1 = num_ent = 0; 00859 ptr = 0; 00860 iMeshP_getEntities( imesh, prtn, part, root, iBase_ALL_TYPES, 00861 iMesh_ALL_TOPOLOGIES, 00862 &ptr, &junk1, &num_ent, &ierr ); CHKERR; 00863 std::vector<iBase_EntityHandle> act_all( ptr, ptr+num_ent ); 00864 free(ptr); 00865 std::sort( part_quads.begin(), part_quads.end() ); 00866 std::sort( part_all.begin(), part_all.end() ); 00867 std::sort( act_quads.begin(), act_quads.end() ); 00868 std::sort( act_all.begin(), act_all.end() ); 00869 ASSERT( part_quads == act_quads ); 00870 ASSERT( part_all == act_all ); 00871 00872 // create an entity set containing half of the quads from the part 00873 std::vector<iBase_EntityHandle> half_quads(part_quads.size()/2); 00874 for (size_t i = 0; i < half_quads.size(); ++i) 00875 half_quads[i] = part_quads[2*i]; 00876 iBase_EntitySetHandle set; 00877 iMesh_createEntSet( imesh, 1, &set, &ierr ); 00878 CHKERR; 00879 iMesh_addEntArrToSet( imesh, &half_quads[0], half_quads.size(), set, &ierr ); 00880 CHKERR; 00881 00882 // check if there exists any quads not in the part that we 00883 // can add to the set 00884 std::vector<iBase_EntityHandle> all_quads, other_quads; 00885 ierr = get_entities( imesh, root, iBase_FACE, iMesh_QUADRILATERAL, all_quads); CHKERR; 00886 std::sort( all_quads.begin(), all_quads.end() ); 00887 std::sort( part_quads.begin(), part_quads.end() ); 00888 std::set_difference( all_quads.begin(), all_quads.end(), 00889 part_quads.begin(), part_quads.end(), 00890 std::back_inserter( other_quads ) ); 00891 iMesh_addEntArrToSet( imesh, &other_quads[0], other_quads.size(), set, &ierr ); 00892 CHKERR; 00893 00894 // compare local counts (using non-root set) 00895 00896 if (test_type) 00897 iMeshP_getNumOfType( imesh, prtn, part, set, iBase_FACE, &count, &ierr ); 00898 else 00899 iMeshP_getNumOfTopo( imesh, prtn, part, set, iMesh_QUADRILATERAL, &count, &ierr ); 00900 CHKERR; 00901 ASSERT( count == (int)half_quads.size() ); 00902 00903 if (test_type) 00904 iMeshP_getNumOfType( imesh, prtn, part, set, iBase_VERTEX, &count, &ierr ); 00905 else 00906 iMeshP_getNumOfTopo( imesh, prtn, part, set, iMesh_POINT, &count, &ierr ); 00907 CHKERR; 00908 ASSERT( count == 0 ); 00909 00910 // compare local contents (using non-root set) 00911 00912 junk1 = 0; num_ent = 0; 00913 ptr = 0; 00914 iMeshP_getEntities( imesh, prtn, part, set, test_type ? iBase_FACE : iBase_ALL_TYPES, 00915 test_type ? iMesh_ALL_TOPOLOGIES : iMesh_QUADRILATERAL, 00916 &ptr, &junk1, &num_ent, &ierr ); CHKERR; 00917 act_quads.resize(num_ent); 00918 std::copy( ptr, ptr + num_ent, act_quads.begin() ); 00919 free(ptr); 00920 std::sort( half_quads.begin(), half_quads.end() ); 00921 std::sort( act_quads.begin(), act_quads.end() ); 00922 ASSERT( act_quads == half_quads ); 00923 00924 return iBase_SUCCESS; 00925 } 00926 00927 00928 00937 int test_get_by_type( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 00938 { 00939 int ierr; 00940 ierr = test_get_by_type_topo_all( imesh, prtn, true, map.num_parts() ); 00941 PCHECK; 00942 ierr = test_get_by_type_topo_local( imesh, prtn, true ); 00943 PCHECK; 00944 return 0; 00945 } 00946 00955 int test_get_by_topo( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 00956 { 00957 int ierr; 00958 ierr = test_get_by_type_topo_all( imesh, prtn, false, map.num_parts() ); 00959 PCHECK; 00960 ierr = test_get_by_type_topo_local( imesh, prtn, false ); 00961 PCHECK; 00962 return 0; 00963 } 00964 00965 00974 int test_part_id_handle( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 00975 { 00976 // get local part ids 00977 int rank, ierr; 00978 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00979 std::vector<iMeshP_Part> ids; 00980 map.part_id_from_rank( rank, ids ); 00981 00982 // check single-part functions and build list of part handles 00983 std::vector<iMeshP_PartHandle> handles( ids.size() ); 00984 size_t i; 00985 for (i = 0; i < ids.size(); ++i) { 00986 iMeshP_getPartHandleFromPartId( imesh, prtn, ids[i], &handles[i], &ierr ); 00987 CHKERR; 00988 iMeshP_Part id; 00989 iMeshP_getPartIdFromPartHandle( imesh, prtn, handles[i], &id, &ierr ); 00990 CHKERR; 00991 if (id != ids[i]) 00992 break; 00993 } 00994 ASSERT( i == ids.size() ); 00995 00996 // test iMeshP_getPartIdsFromPartHandlesArr 00997 std::vector<iMeshP_Part> ids2( ids.size() ); 00998 int junk1 = ids.size(), junk2 = 0; 00999 iMeshP_Part* ptr = &ids2[0]; 01000 iMeshP_getPartIdsFromPartHandlesArr( imesh, prtn, &handles[0], handles.size(), 01001 &ptr, &junk1, &junk2, &ierr ); 01002 PCHECK; 01003 ASSERT( ptr == &ids2[0] ); 01004 ASSERT( junk2 == (int)ids2.size() ); 01005 ASSERT( ids == ids2 ); 01006 01007 // test iMeshP_getPartHandlesFromPartsIdsArr 01008 std::vector<iMeshP_PartHandle> handles2(handles.size()); 01009 junk1 = handles.size(); 01010 junk2 = 0; 01011 iMeshP_PartHandle* ptr2 = &handles2[0]; 01012 iMeshP_getPartHandlesFromPartsIdsArr( imesh, prtn, &ids[0], ids.size(), 01013 &ptr2, &junk1, &junk2, &ierr ); 01014 PCHECK; 01015 ASSERT( ptr2 == &handles2[0] ); 01016 ASSERT( junk2 == (int)handles2.size() ); 01017 ASSERT( handles == handles2 ); 01018 01019 return 0; 01020 } 01021 01028 int test_part_rank( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 01029 { 01030 int ierr = 0, rank; 01031 std::vector<iMeshP_Part> invalid, failed; 01032 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01033 01034 // test iMeshP_getRankOfPart 01035 for (size_t i = 0; i < map.get_parts().size(); ++i) { 01036 int pr; 01037 iMeshP_getRankOfPart( imesh, prtn, map.get_parts()[i], &pr, &ierr ); 01038 if (iBase_SUCCESS != ierr) 01039 failed.push_back( map.get_parts()[i] ); 01040 else if (pr != map.get_ranks()[i]) 01041 invalid.push_back( map.get_parts()[i] ); 01042 } 01043 if (!failed.empty()) { 01044 std::cerr << "Processor " << rank << ": iMeshP_getRankOfPart failed for " << failed.size() << " parts." << std::endl; 01045 ierr = iBase_FAILURE; 01046 } 01047 if (!invalid.empty()) { 01048 std::cerr << "Processor " << rank << ": iMeshP_getRankOfPart was incorrect for " << invalid.size() << " parts." << std::endl; 01049 ierr = iBase_FAILURE; 01050 } 01051 PCHECK; 01052 01053 // test iMeshP_getRankOfPartArr 01054 std::vector<int> ranks( map.get_parts().size() ); 01055 int junk1 = ranks.size(), junk2, *ptr = &ranks[0]; 01056 iMeshP_getRankOfPartArr( imesh, prtn, &map.get_parts()[0], map.get_parts().size(), 01057 &ptr, &junk1, &junk2, &ierr ); 01058 PCHECK; 01059 assert( ptr == &ranks[0] ); 01060 assert( junk1 == (int)ranks.size() ); 01061 ASSERT( junk2 == (int)ranks.size() ); 01062 for (size_t i = 0; i < map.get_parts().size(); ++i) { 01063 if (ranks[i] != map.get_ranks()[i]) 01064 invalid.push_back( map.get_parts()[i] ); 01065 } 01066 if (!invalid.empty()) { 01067 std::cerr << "Processor " << rank << ": iMeshP_getRankOfPartArr was incorrect for " << invalid.size() << " parts." << std::endl; 01068 ierr = iBase_FAILURE; 01069 } 01070 PCHECK; 01071 01072 return 0; 01073 } 01074 01075 01076 // see create_mesh(..) 01077 static void get_part_neighbors( int logical_part_id, 01078 int num_parts, 01079 int neighbors[5], 01080 int& num_neighbors ) 01081 { 01082 num_neighbors = 0; 01083 if (logical_part_id + 1 < num_parts) 01084 neighbors[num_neighbors++] = logical_part_id + 1; 01085 if (logical_part_id + 2 < num_parts) 01086 neighbors[num_neighbors++] = logical_part_id + 2; 01087 if (logical_part_id % 2) { 01088 neighbors[num_neighbors++] = logical_part_id - 1; 01089 if (logical_part_id > 2) { 01090 neighbors[num_neighbors++] = logical_part_id - 3; 01091 neighbors[num_neighbors++] = logical_part_id - 2; 01092 } 01093 } 01094 else { 01095 if (logical_part_id + 3 < num_parts) 01096 neighbors[num_neighbors++] = logical_part_id + 3; 01097 if (logical_part_id > 1) { 01098 neighbors[num_neighbors++] = logical_part_id - 1; 01099 neighbors[num_neighbors++] = logical_part_id - 2; 01100 } 01101 } 01102 } 01103 01112 int test_get_neighbors( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 01113 { 01114 int ierr, rank; 01115 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01116 01117 std::vector<iMeshP_Part> local_parts; 01118 map.part_id_from_rank( rank, local_parts ); 01119 01120 // get handles for local parts 01121 std::vector<iMeshP_PartHandle> handles(local_parts.size()); 01122 iMeshP_PartHandle* ptr = &handles[0]; 01123 int junk1 = handles.size(), junk2 = 0; 01124 iMeshP_getPartHandlesFromPartsIdsArr( imesh, prtn, &local_parts[0], local_parts.size(), 01125 &ptr, &junk1, &junk2, &ierr ); 01126 PCHECK; 01127 assert( ptr == &handles[0] ); 01128 assert( junk2 == (int)handles.size() ); 01129 01130 // get logical ids for local parts 01131 std::vector<int> logical_ids; 01132 map.local_id_from_rank( rank, logical_ids ); 01133 01134 // get neighbors for each local part 01135 std::vector< std::vector<iMeshP_Part> > neighbors( logical_ids.size() ); 01136 for (size_t i = 0;i < logical_ids.size(); ++i) { 01137 int logical_neighbors[5], num_neighbors; 01138 get_part_neighbors( logical_ids[i], map.num_parts(), logical_neighbors, num_neighbors ); 01139 neighbors[i].resize( num_neighbors ); 01140 for (int j = 0; j < num_neighbors; ++j) 01141 neighbors[i][j] = map.part_id_from_local_id( logical_neighbors[j] ); 01142 std::sort( neighbors[i].begin(), neighbors[i].end() ); 01143 } 01144 01145 // test iMeshP_getNumPartNbors 01146 std::vector< iMeshP_Part > invalid, failed; 01147 for (size_t i = 0; i < local_parts.size(); ++i) { 01148 int count; 01149 iMeshP_getNumPartNbors( imesh, prtn, handles[i], iBase_VERTEX, &count, &ierr ); 01150 if (ierr) 01151 failed.push_back( local_parts[i] ); 01152 else if (count != (int)neighbors[i].size()) 01153 invalid.push_back( local_parts[i] ); 01154 } 01155 if (!failed.empty()) { 01156 std::cerr << "Processor " << rank << ": iMeshP_getNumPartNbors failed for " << failed.size() << " parts." << std::endl; 01157 ierr = iBase_FAILURE; PCHECK; 01158 } 01159 if (!invalid.empty()) { 01160 std::cerr << "Processor " << rank << ": iMeshP_getNumPartNbors was incorrect for " << invalid.size() << " parts." << std::endl; 01161 ierr = iBase_FAILURE; PCHECK; 01162 } 01163 01164 // test iMeshP_getPartNbors 01165 ierr = 0; 01166 for (size_t i = 0; i < local_parts.size(); ++i) { 01167 int count, junk = 0, another_count; 01168 iMeshP_Part* list = 0; 01169 iMeshP_getPartNbors( imesh, prtn, handles[i], iBase_VERTEX, &another_count, &list, &junk, &count, &ierr ); 01170 assert( count == another_count ); 01171 if (ierr) 01172 failed.push_back( local_parts[i] ); 01173 else { 01174 std::sort( list, list+count ); 01175 std::vector<iMeshP_Part> cpy( list, list+count ); 01176 if (cpy != neighbors[i]) 01177 invalid.push_back( local_parts[i] ); 01178 free(list); 01179 } 01180 } 01181 if (!failed.empty()) { 01182 std::cerr << "Processor " << rank << ": iMeshP_getPartNbors failed for " << failed.size() << " parts." << std::endl; 01183 ierr = iBase_FAILURE; 01184 } 01185 if (!invalid.empty()) { 01186 std::cerr << "Processor " << rank << ": iMeshP_getPartNbors was incorrect for " << invalid.size() << " parts." << std::endl; 01187 ierr = iBase_FAILURE; 01188 } 01189 PCHECK; 01190 01191 // test iMeshP_getNumPartNborsArr 01192 std::vector<int> count_vect( handles.size() ); 01193 int* count_arr = &count_vect[0]; 01194 junk1 = handles.size(); 01195 iMeshP_getNumPartNborsArr( imesh, prtn, &handles[0], handles.size(), iBase_VERTEX, 01196 &count_arr, &junk1, &junk2, &ierr ); 01197 PCHECK; 01198 assert( count_arr = &count_vect[0] ); 01199 assert( junk2 == (int)handles.size() ); 01200 for (size_t i = 0; i < local_parts.size(); ++i) { 01201 if (count_arr[i] != (int)neighbors[i].size()) 01202 invalid.push_back( local_parts[i] ); 01203 } 01204 if (!invalid.empty()) { 01205 std::cerr << "Processor " << rank << ": iMeshP_getNumPartNborsArr was incorrect for " << invalid.size() << " parts." << std::endl; 01206 ierr = iBase_FAILURE; 01207 } 01208 PCHECK; 01209 01210 // test iMeshP_getPartNborsArr 01211 iMeshP_Part* nbor_arr = 0; 01212 junk1 = handles.size(), junk2 = 0; 01213 int junk3 = 0, nbor_size; 01214 iMeshP_getPartNborsArr( imesh, prtn, &handles[0], handles.size(), iBase_VERTEX, 01215 &count_arr, &junk1, &junk2, 01216 &nbor_arr, &junk3, &nbor_size, 01217 &ierr ); 01218 PCHECK; 01219 assert( count_arr = &count_vect[0] ); 01220 assert( junk2 == (int)handles.size() ); 01221 std::vector<iMeshP_Part> all_nbors( nbor_arr, nbor_arr + nbor_size ); 01222 free( nbor_arr ); 01223 std::vector<iMeshP_Part>::iterator j = all_nbors.begin(); 01224 bool bad_length = false; 01225 for (size_t i = 0; i < local_parts.size(); ++i) { 01226 if (all_nbors.end() - j > count_arr[i]) { 01227 bad_length = true; 01228 break; 01229 } 01230 if (count_arr[i] != (int)neighbors[i].size()) { 01231 invalid.push_back( local_parts[i] ); 01232 } 01233 else { 01234 std::vector<iMeshP_Part>::iterator e = j + count_arr[i]; 01235 std::sort( j, e ); 01236 if (!std::equal( j, e, neighbors[i].begin() )) 01237 invalid.push_back( local_parts[i] ); 01238 } 01239 } 01240 if (bad_length) { 01241 std::cerr << "Processor " << rank << ": iMeshP_getPartNborsArr had inconsistent result array lengths." << std::endl; 01242 ierr = iBase_FAILURE; 01243 } 01244 if (!invalid.empty()) { 01245 std::cerr << "Processor " << rank << ": iMeshP_getPartNborsArr was incorrect for " << invalid.size() << " parts." << std::endl; 01246 ierr = iBase_FAILURE; 01247 } 01248 PCHECK; 01249 01250 return 0; 01251 } 01252 01253 // Determine the expected vertices on the interface between two parts. 01254 // Returns no vertices for non-adjacient parts and fails if both parts 01255 // are the same. 01256 // See create_mesh(..) for the assumed mesh. 01257 static int interface_verts( iMesh_Instance imesh, 01258 iMeshP_PartitionHandle prtn, 01259 iMeshP_PartHandle local_part, 01260 iMeshP_Part other_part, 01261 const PartMap& map, 01262 std::vector<iBase_EntityHandle> &vtx_handles ) 01263 { 01264 int ierr, rank; 01265 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01266 01267 iMeshP_Part local_id; 01268 iMeshP_getPartIdFromPartHandle( imesh, prtn, local_part, &local_id, &ierr ); 01269 CHKERR; 01270 01271 const int local_logical = map.local_id_from_part_id( local_id ); 01272 const int other_logical = map.local_id_from_part_id( other_part ); 01273 01274 // get grid of local vertices 01275 01276 iBase_EntityHandle verts[3][3]; 01277 const double xbase = (local_id / 2) * 2; 01278 const double ybase = (local_id % 2) * 2; 01279 01280 // get quads in partition 01281 iBase_EntityHandle quads[4], *ptr = quads; 01282 int junk1 = 4, junk2; 01283 iMesh_getEntities( imesh, local_part, iBase_FACE, iMesh_QUADRILATERAL, &ptr, &junk1, &junk2, &ierr ); 01284 CHKERR; 01285 assert( ptr == quads ); 01286 assert( junk1 == 4 ); 01287 assert( junk2 == 4 ); 01288 01289 // get vertices in quads 01290 iBase_EntityHandle conn[16]; 01291 int offsets[5], *off_ptr = offsets, junk3 = 5, junk4; 01292 ptr = conn; 01293 junk1 = 16; 01294 iMesh_getEntArrAdj( imesh, quads, 4, iBase_VERTEX, 01295 &ptr, &junk1, &junk2, 01296 &off_ptr, &junk3, &junk4, 01297 &ierr ); 01298 CHKERR; 01299 assert( ptr == conn ); 01300 assert( junk1 == 16 ); 01301 assert( junk2 == 16 ); 01302 assert( off_ptr == offsets ); 01303 assert( junk3 == 5 ); 01304 assert( junk4 == 5 ); 01305 01306 // make unique vertex list 01307 std::sort( conn, conn + 16 ); 01308 const int num_vtx = std::unique( conn, conn+16 ) - conn; 01309 assert(9 == num_vtx); 01310 01311 // get vertex coords 01312 std::vector<double> coords(27); 01313 ierr = get_coords( imesh, conn, 9, &coords[0] ); 01314 CHKERR; 01315 01316 // use vertex coords to determine logical position 01317 for (int i = 0; i < num_vtx; ++i) { 01318 int x = (int)round(coords[3*i ] - xbase); 01319 int y = (int)round(coords[3*i+1] - ybase); 01320 if (x < 0 || x > 2 || y < 0 || y > 2) { 01321 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01322 << " Invalid vertex coordinate: (" << coords[3*i] << ", " << coords[3*i+1] 01323 << ", " << coords[3*i+2] << ")" << std::endl 01324 << " For logical partition " << local_id << std::endl; 01325 return iBase_FAILURE; 01326 } 01327 verts[x][y] = conn[i]; 01328 } 01329 01330 if (local_logical % 2) { 01331 switch (other_logical - local_logical) { 01332 case 0: 01333 return iBase_FAILURE; 01334 case 1: // upper right 01335 vtx_handles.resize(1); 01336 vtx_handles[0] = verts[2][0]; 01337 break; 01338 case 2: // right 01339 vtx_handles.resize(3); 01340 std::copy( verts[2], verts[2]+3, vtx_handles.begin() ); 01341 break; 01342 case -1: // above 01343 vtx_handles.resize(3); 01344 vtx_handles[0] = verts[0][0]; 01345 vtx_handles[1] = verts[1][0]; 01346 vtx_handles[2] = verts[2][0]; 01347 break; 01348 case -2: // left 01349 vtx_handles.resize(3); 01350 std::copy( verts[0], verts[0]+3, vtx_handles.begin() ); 01351 break; 01352 case -3: // upper left 01353 vtx_handles.resize(1); 01354 vtx_handles[0] = verts[0][0]; 01355 break; 01356 default: 01357 vtx_handles.clear(); 01358 break; 01359 } 01360 } 01361 else { 01362 switch (other_logical - local_logical) { 01363 case 0: 01364 return iBase_FAILURE; 01365 case 1: // below 01366 vtx_handles.resize(3); 01367 vtx_handles[0] = verts[0][2]; 01368 vtx_handles[1] = verts[1][2]; 01369 vtx_handles[2] = verts[2][2]; 01370 break; 01371 case 2: // right 01372 vtx_handles.resize(3); 01373 std::copy( verts[2], verts[2]+3, vtx_handles.begin() ); 01374 break; 01375 case 3: // lower right 01376 vtx_handles.resize(1); 01377 vtx_handles[0] = verts[2][2]; 01378 break; 01379 case -1: // lower left 01380 vtx_handles.resize(1); 01381 vtx_handles[0] = verts[0][2]; 01382 break; 01383 case -2: // left 01384 vtx_handles.resize(3); 01385 std::copy( verts[0], verts[0]+3, vtx_handles.begin() ); 01386 break; 01387 default: 01388 vtx_handles.clear(); 01389 break; 01390 } 01391 } 01392 01393 return iBase_SUCCESS; 01394 } 01395 01396 01397 01398 01399 01406 int test_get_part_boundary( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 01407 { 01408 int ierr, rank; 01409 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01410 01411 // get local part handles and part ids, and global part id list 01412 std::vector<iMeshP_PartHandle> local_handles; 01413 std::vector<iMeshP_Part> local_ids; 01414 std::vector<iMeshP_Part> all_parts = map.get_parts(); 01415 std::map< iMeshP_PartHandle, std::vector<iBase_EntityHandle> > part_bdry; 01416 ierr = get_local_parts( imesh, prtn, local_handles, &local_ids ); 01417 CHKERR; 01418 01419 // for each combination of local part with any other part, 01420 // check for valid function values. 01421 std::vector< std::pair<iMeshP_Part,iMeshP_Part> > num_failed, num_error, list_failed, list_error, error; 01422 for (size_t i = 0; i < local_handles.size(); ++i) { 01423 iMeshP_PartHandle local_handle = local_handles[i]; 01424 iMeshP_Part local_id = local_ids[i]; 01425 for (std::vector<iMeshP_Part>::iterator j = all_parts.begin(); j != all_parts.end(); ++j) { 01426 iMeshP_Part other_id = *j; 01427 if (other_id == local_id) 01428 continue; 01429 01430 std::pair<iMeshP_Part,iMeshP_Part> part_pair; 01431 part_pair.first = local_id; 01432 part_pair.second = other_id; 01433 01434 // get expected values 01435 std::vector<iBase_EntityHandle> shared_verts; 01436 ierr = interface_verts( imesh, prtn, local_handle, other_id, map, shared_verts ); 01437 if (ierr != iBase_SUCCESS) { 01438 error.push_back( part_pair ); 01439 continue; 01440 } 01441 std::sort( shared_verts.begin(), shared_verts.end() ); 01442 01443 // test iMeshP_getNumPartBdryEnts 01444 int count; 01445 iMeshP_getNumPartBdryEnts( imesh, prtn, local_handle, iBase_VERTEX, iMesh_POINT, 01446 other_id, &count, &ierr ); 01447 if (iBase_SUCCESS != ierr) 01448 num_error.push_back( part_pair ); 01449 else if (count != (int)shared_verts.size()) 01450 num_failed.push_back( part_pair ); 01451 01452 // test iMeshP_getPartBdryEnts 01453 iBase_EntityHandle* ptr = 0; 01454 int junk = 0; 01455 iMeshP_getPartBdryEnts( imesh, prtn, local_handle, iBase_VERTEX, iMesh_POINT, other_id, 01456 &ptr, &junk, &count, &ierr ); 01457 if (iBase_SUCCESS != ierr) 01458 list_error.push_back( part_pair ); 01459 else { 01460 std::copy( ptr, ptr + count, std::back_inserter( part_bdry[local_handles[i]] ) ); 01461 std::sort( ptr, ptr + count ); 01462 if ((int)shared_verts.size() != count || 01463 !std::equal( shared_verts.begin(), shared_verts.end(), ptr )) 01464 list_failed.push_back( part_pair ); 01465 free(ptr); 01466 } 01467 } 01468 } 01469 01470 if (!error.empty()) { 01471 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01472 << " Internal error for " << error.size() << " part pairs." << std::endl; 01473 ierr = iBase_FAILURE; 01474 } 01475 if (!num_error.empty()) { 01476 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01477 << " iMeshP_getNumPartBdryEnts return error for " << num_error.size() << " part pairs." << std::endl; 01478 ierr = iBase_FAILURE; 01479 } 01480 if (!list_error.empty()) { 01481 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01482 << " iMeshP_getPartBdryEnts return error for " << list_error.size() << " part pairs." << std::endl; 01483 ierr = iBase_FAILURE; 01484 } 01485 if (!num_failed.empty()) { 01486 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01487 << " iMeshP_getNumPartBdryEnts return incorrect results for " << num_failed.size() << " part pairs." << std::endl; 01488 ierr = iBase_FAILURE; 01489 } 01490 if (!list_failed.empty()) { 01491 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01492 << " iMeshP_getPartBdryEnts return incorrect results for " << list_failed.size() << " part pairs." << std::endl; 01493 ierr = iBase_FAILURE; 01494 } 01495 01496 if (iBase_SUCCESS != ierr) 01497 return ierr; 01498 01499 01500 // test with iMeshP_ALL_PARTS 01501 for (size_t i = 0; i < local_handles.size(); ++i) { 01502 std::vector<iBase_EntityHandle>& exp_bdry = part_bdry[local_handles[i]]; 01503 std::sort( exp_bdry.begin(), exp_bdry.end() ); 01504 exp_bdry.erase( std::unique( exp_bdry.begin(), exp_bdry.end() ), exp_bdry.end() ); 01505 std::pair<iMeshP_Part,iMeshP_Part> part_pair; 01506 part_pair.first = local_ids[i]; 01507 part_pair.second = iMeshP_ALL_PARTS; 01508 01509 int num = 0; 01510 iMeshP_getNumPartBdryEnts( imesh, prtn, local_handles[i], iBase_VERTEX, iMesh_POINT, 01511 iMeshP_ALL_PARTS, &num, &ierr ); 01512 if (ierr) 01513 num_error.push_back( part_pair ); 01514 else if (num != (int)exp_bdry.size()) 01515 num_failed.push_back( part_pair ); 01516 01517 iBase_EntityHandle* bdry = 0; 01518 int junk = num = 0; 01519 iMeshP_getPartBdryEnts( imesh, prtn, local_handles[i], iBase_VERTEX, iMesh_POINT, 01520 iMeshP_ALL_PARTS, &bdry, &junk, &num, &ierr ); 01521 if (ierr) 01522 list_error.push_back( part_pair ); 01523 else { 01524 std::sort( bdry, bdry+num ); 01525 if (num != (int)exp_bdry.size() || !std::equal( bdry, bdry+num, exp_bdry.begin() )) 01526 list_failed.push_back( part_pair ); 01527 free(bdry); 01528 } 01529 } 01530 if (!num_error.empty()) { 01531 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01532 << " iMeshP_getNumPartBdryEnts return error for " << num_error.size() << " part pairs." << std::endl; 01533 ierr = iBase_FAILURE; 01534 } 01535 if (!list_error.empty()) { 01536 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01537 << " iMeshP_getPartBdryEnts return error for " << list_error.size() << " part pairs." << std::endl; 01538 ierr = iBase_FAILURE; 01539 } 01540 if (!num_failed.empty()) { 01541 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01542 << " iMeshP_getNumPartBdryEnts return incorrect results for " << num_failed.size() << " part pairs." << std::endl; 01543 ierr = iBase_FAILURE; 01544 } 01545 if (!list_failed.empty()) { 01546 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01547 << " iMeshP_getPartBdryEnts return incorrect results for " << list_failed.size() << " part pairs." << std::endl; 01548 ierr = iBase_FAILURE; 01549 } 01550 01551 return ierr; 01552 } 01553 01560 int test_part_boundary_iter( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 01561 { 01562 int ierr, rank, has_data; 01563 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01564 01565 // get local part handles and part ids, and global part id list 01566 std::vector<iMeshP_PartHandle> local_handles; 01567 std::vector<iMeshP_Part> local_ids; 01568 std::vector<iMeshP_Part> all_parts = map.get_parts(); 01569 ierr = get_local_parts( imesh, prtn, local_handles, &local_ids ); 01570 CHKERR; 01571 01572 std::vector< std::pair<iMeshP_Part,iMeshP_Part> > single_failed, single_error, 01573 single_step_error, array_failed, 01574 array_error, array_step_error; 01575 for (size_t i = 0; i < local_handles.size(); ++i) { 01576 iMeshP_PartHandle local_handle = local_handles[i]; 01577 iMeshP_Part local_id = local_ids[i]; 01578 for (std::vector<iMeshP_Part>::iterator j = all_parts.begin(); j != all_parts.end(); ++j) { 01579 iMeshP_Part other_id = *j; 01580 if (other_id == local_id) 01581 continue; 01582 01583 std::pair<iMeshP_Part,iMeshP_Part> part_pair; 01584 part_pair.first = local_id; 01585 part_pair.second = other_id; 01586 01587 // get expected values 01588 std::vector<iBase_EntityHandle> shared_verts; 01589 ierr = interface_verts( imesh, prtn, local_handle, other_id, map, shared_verts ); 01590 if (ierr != iBase_SUCCESS || 0 == shared_verts.size()) 01591 continue; 01592 std::sort( shared_verts.begin(), shared_verts.end() ); 01593 01594 // test single entity iterator 01595 iBase_EntityIterator siter; 01596 iMeshP_initPartBdryEntIter( imesh, prtn, local_handle, iBase_VERTEX, iMesh_POINT, 01597 other_id, &siter, &ierr ); 01598 if (ierr != iBase_SUCCESS) { 01599 single_error.push_back( part_pair ); 01600 } 01601 else { 01602 std::vector<iBase_EntityHandle> results; 01603 for (;;) { 01604 iBase_EntityHandle handle; 01605 iMesh_getNextEntIter( imesh, siter, &handle, &has_data, &ierr ); 01606 if (ierr != iBase_SUCCESS) { 01607 single_step_error.push_back( part_pair ); 01608 break; 01609 } 01610 if (!has_data) 01611 break; 01612 results.push_back(handle); 01613 } 01614 01615 std::sort( results.begin(), results.end() ); 01616 if (results.size() != shared_verts.size() || 01617 !std::equal( results.begin(), results.end(), shared_verts.begin())) 01618 single_failed.push_back( part_pair ); 01619 } 01620 iMesh_endEntIter( imesh, siter, &ierr ); 01621 01622 // test array iterator 01623 iBase_EntityArrIterator aiter; 01624 iMeshP_initPartBdryEntArrIter( imesh, prtn, local_handle, iBase_VERTEX, iMesh_POINT, 01625 shared_verts.size(), other_id, &aiter, &ierr ); 01626 if (ierr != iBase_SUCCESS) { 01627 array_error.push_back( part_pair ); 01628 continue; 01629 } 01630 iBase_EntityHandle results[5], *ptr = results; 01631 int junk = 5, count; 01632 iMesh_getNextEntArrIter( imesh, aiter, &ptr, &junk, &count, &has_data, &ierr ); 01633 if (ierr != iBase_SUCCESS || !has_data) { 01634 array_step_error.push_back( part_pair ); 01635 continue; 01636 } 01637 assert(count <= 5); 01638 assert(ptr == results); 01639 std::sort(ptr, ptr + count); 01640 if (count != (int)shared_verts.size() || 01641 !std::equal( shared_verts.begin(), shared_verts.end(), results )) 01642 array_failed.push_back( part_pair ); 01643 } 01644 } 01645 01646 if (!single_error.empty()) { 01647 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01648 << " iMeshP_initPartBdryEntIter return error for " << single_error.size() << " part pairs." << std::endl; 01649 ierr = iBase_FAILURE; 01650 } 01651 if (!single_step_error.empty()) { 01652 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01653 << " iMesh_getNextEntIter return error for " << single_step_error.size() << " part pairs." << std::endl; 01654 ierr = iBase_FAILURE; 01655 } 01656 if (!single_failed.empty()) { 01657 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01658 << " iMeshP_initPartBdryEntIter iterator iterated over invalid entities for " 01659 << single_failed.size() << " part pairs." << std::endl; 01660 ierr = iBase_FAILURE; 01661 } 01662 01663 if (!array_error.empty()) { 01664 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01665 << " iMeshP_initPartBdryEntArrIter return error for " << array_error.size() << " part pairs." << std::endl; 01666 ierr = iBase_FAILURE; 01667 } 01668 if (!array_step_error.empty()) { 01669 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01670 << " iMesh_getNextEntArrIter return error for " << array_step_error.size() << " part pairs." << std::endl; 01671 ierr = iBase_FAILURE; 01672 } 01673 if (!array_failed.empty()) { 01674 std::cerr << "Processor " << rank << ": Error at " __FILE__ ":" << __LINE__ << std::endl 01675 << " iMeshP_initPartBdryEntArrIter iterator iterated over invalid entities for " 01676 << array_failed.size() << " part pairs." << std::endl; 01677 ierr = iBase_FAILURE; 01678 } 01679 01680 return ierr; 01681 } 01682 01688 int test_get_adjacencies( iMesh_Instance /* imesh */, iMeshP_PartitionHandle /* prtn */, const PartMap& ) 01689 { 01690 return iBase_SUCCESS; 01691 } 01692 01699 int test_entity_iterator( iMesh_Instance /*imesh */, iMeshP_PartitionHandle /*prtn*/, const PartMap& ) 01700 { 01701 return iBase_SUCCESS; 01702 } 01703 01712 int test_entity_owner( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& /* map */ ) 01713 { 01714 int ierr, rank, size; 01715 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01716 MPI_Comm_size( MPI_COMM_WORLD, &size ); 01717 01718 // get local part handles and part ids 01719 std::vector<iMeshP_PartHandle> local_handles; 01720 std::vector<iMeshP_Part> local_ids; 01721 ierr = get_local_parts( imesh, prtn, local_handles, &local_ids ); 01722 PCHECK; 01723 01724 // test iMeshP_getEntOwnerPart for quads in each part 01725 std::vector<iBase_EntityHandle> all_quads; 01726 std::vector<iMeshP_Part> quad_owners; 01727 int invalid_count = 0; 01728 for (size_t i = 0; i < local_handles.size(); ++i) { 01729 std::vector<iBase_EntityHandle> quads; 01730 ierr = get_entities( imesh, local_handles[0], iBase_FACE, iMesh_QUADRILATERAL, quads ); 01731 if (ierr) 01732 break; 01733 01734 for (size_t j = 0; j < quads.size(); ++j) { 01735 all_quads.push_back( quads[j] ); 01736 quad_owners.push_back( local_ids[i] ); 01737 iMeshP_Part owner; 01738 iMeshP_getEntOwnerPart( imesh, prtn, quads[j], &owner, &ierr ); 01739 if (iBase_SUCCESS != ierr) 01740 break; 01741 01742 if (owner != local_ids[i]) 01743 ++invalid_count; 01744 } 01745 if (iBase_SUCCESS != ierr) 01746 break; 01747 } 01748 PCHECK; 01749 ASSERT(0 == invalid_count); 01750 01751 // test iMeshP_getEntOwnerPartArr for quads in each part 01752 invalid_count = 0; 01753 for (size_t i = 0; i < local_handles.size(); ++i) { 01754 std::vector<iBase_EntityHandle> quads; 01755 ierr = get_entities( imesh, local_handles[0], iBase_FACE, iMesh_QUADRILATERAL, quads ); 01756 if (ierr) 01757 break; 01758 01759 std::vector<iMeshP_Part> owners(quads.size()), expected(quads.size(),local_ids[i]); 01760 int junk = owners.size(), count; 01761 iMeshP_Part* ptr = &owners[0]; 01762 iMeshP_getEntOwnerPartArr( imesh, prtn, &quads[0], quads.size(), 01763 &ptr, &junk, &count, &ierr ); 01764 if (ierr) 01765 break; 01766 assert( ptr == &owners[0] ); 01767 assert( junk == (int)owners.size() ); 01768 assert( count == (int)quads.size() ); 01769 if (owners != expected) 01770 ++invalid_count; 01771 } 01772 PCHECK; 01773 ASSERT(0 == invalid_count); 01774 01775 // get all vertices 01776 iBase_EntityHandle* vtx_arr = 0; 01777 int junk1 = 0, num_vtx; 01778 int *junk2 = 0, junk3 = 0, junk4; 01779 iMesh_getEntArrAdj( imesh, &all_quads[0], all_quads.size(), iBase_VERTEX, 01780 &vtx_arr, &junk1, &num_vtx, &junk2, &junk3, &junk4, &ierr ); 01781 PCHECK; 01782 free(junk2); 01783 std::sort( vtx_arr, vtx_arr + num_vtx ); 01784 num_vtx = std::unique( vtx_arr, vtx_arr + num_vtx ) - vtx_arr; 01785 std::vector<iBase_EntityHandle> all_verts( vtx_arr, vtx_arr + num_vtx ); 01786 free(vtx_arr); 01787 01788 // check consistency between iMeshP_getEntOwnerPart and iMeshP_getEntOwnerPartArr 01789 // for all vertices 01790 std::vector<iMeshP_Part> vert_owners( all_verts.size() ); 01791 junk1 = vert_owners.size(); 01792 iMeshP_Part* junk5 = &vert_owners[0]; 01793 iMeshP_getEntOwnerPartArr( imesh, prtn, &all_verts[0], all_verts.size(), 01794 &junk5, &junk1, &junk3, &ierr ); 01795 PCHECK; 01796 assert( junk5 == &vert_owners[0] ); 01797 assert( junk1 == (int)vert_owners.size() ); 01798 assert( junk3 == (int)all_verts.size() ); 01799 01800 invalid_count = 0; 01801 for (size_t i = 0; i < all_verts.size(); ++i) { 01802 iMeshP_Part owner; 01803 iMeshP_getEntOwnerPart( imesh, prtn, all_verts[i], &owner, &ierr ); 01804 if (iBase_SUCCESS != ierr || owner != vert_owners[i]) 01805 ++invalid_count; 01806 } 01807 ASSERT(0 == invalid_count); 01808 01809 // get lists for all entities 01810 std::vector<iBase_EntityHandle> all_entities(all_verts); 01811 std::copy( all_quads.begin(), all_quads.end(), std::back_inserter(all_entities) ); 01812 std::vector<iMeshP_Part> all_owners( vert_owners ); 01813 std::copy( quad_owners.begin(), quad_owners.end(), std::back_inserter(all_owners) ); 01814 01815 // check consistency of iMeshP_isEntOwner for all entities 01816 invalid_count = 0; 01817 ierr = iBase_SUCCESS; 01818 for (size_t i = 0; i < local_handles.size(); ++i) { 01819 for (size_t j = 0; ierr == iBase_SUCCESS && j < all_entities.size(); ++j) { 01820 int is_owner; 01821 iMeshP_isEntOwner( imesh, prtn, local_handles[i], all_entities[j], &is_owner, &ierr ); 01822 if (ierr != iBase_SUCCESS) 01823 break; 01824 if (!is_owner == (local_ids[i] == all_owners[j])) 01825 ++invalid_count; 01826 } 01827 } 01828 PCHECK; 01829 ASSERT(0 == invalid_count); 01830 01831 // check consistency of iMeshP_isEntOwnerArr for all entities 01832 for (size_t i = 0; i < local_handles.size(); ++i) { 01833 std::vector<int> is_owner_list(all_entities.size()); 01834 junk1 = is_owner_list.size(); 01835 int* junk6 = &is_owner_list[0]; 01836 iMeshP_isEntOwnerArr( imesh, prtn, local_handles[i], &all_entities[0], all_entities.size(), 01837 &junk6, &junk1, &junk3, &ierr ); 01838 if (iBase_SUCCESS != ierr) 01839 break; 01840 assert( junk6 == &is_owner_list[0] ); 01841 assert( junk1 == (int)is_owner_list.size() ); 01842 assert( junk3 == (int)all_entities.size() ); 01843 invalid_count = 0; 01844 for (size_t j = 0; j < all_entities.size(); ++j) { 01845 if (!(is_owner_list[j]) == (local_ids[0] == all_owners[j])) 01846 ++invalid_count; 01847 } 01848 } 01849 PCHECK; 01850 ASSERT(0 == invalid_count); 01851 01852 01853 // check globally consistent owners for all vertices 01854 01855 // first communicate total number of vertex entries to be sent to root proc 01856 int local_count = all_verts.size(), global_count = 0; 01857 ierr = MPI_Reduce( &local_count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD ); 01858 CHKERR; 01859 01860 // for each vertex, store { (x << 2) | y, owning part id } 01861 std::vector<int> vtxdata( 2 * all_verts.size() ); 01862 std::vector<double> coords( 3 * all_verts.size() ); 01863 ierr = get_coords( imesh, &all_verts[0], all_verts.size(), &coords[0] ); 01864 CHKERR; 01865 for (size_t i = 0; i < all_verts.size(); ++i) { 01866 int x = (int)round(coords[3*i ]); 01867 int y = (int)round(coords[3*i+1]); 01868 vtxdata[2*i ] = (x << 2) | y; 01869 vtxdata[2*i+1] = vert_owners[i]; 01870 } 01871 01872 // collect all data on root procesor 01873 std::vector<int> all_data( 2*global_count ); 01874 std::vector<int> displ(size), counts(size); 01875 if (1 == size) { 01876 std::copy(vtxdata.begin(), vtxdata.end(), all_data.begin()); 01877 counts[0] = vtxdata.size(); 01878 displ[0] = 0; 01879 } 01880 else { 01881 ierr = MPI_Gatherv( &vtxdata[0], vtxdata.size(), MPI_INT, 01882 &all_data[0], &counts[0], &displ[0], MPI_INT, 01883 0, MPI_COMM_WORLD ); 01884 CHKERR; 01885 } 01886 if (rank == 0) { 01887 // map from vertex tag to indices into data 01888 std::multimap<int,int> data_map; 01889 for (int i = 0; i < global_count; ++i) { 01890 std::pair<int,int> p; 01891 p.first = all_data[2*i]; 01892 p.second = i; 01893 data_map.insert( p ); 01894 } 01895 01896 // check consistent data for each vtx 01897 std::multimap<int,int>::const_iterator a, b; 01898 for (a = data_map.begin(); a != data_map.end(); a = b) { 01899 for (b = a; b != data_map.end() && a->first == b->first; ++b) { 01900 int idx1 = a->second; 01901 int idx2 = b->second; 01902 if (all_data[2*idx1+1] == all_data[2*idx2+1]) 01903 continue; 01904 01905 ierr = iBase_FAILURE; 01906 01907 int proc1 = std::lower_bound(displ.begin(), displ.end(),2*idx1) - displ.begin(); 01908 if (displ[proc1] != 2*idx1) 01909 ++proc1; 01910 int proc2 = std::lower_bound(displ.begin(), displ.end(),2*idx2) - displ.begin(); 01911 if (displ[proc2] != 2*idx2) 01912 ++proc2; 01913 01914 std::cerr << "Error at " __FILE__ ":" << __LINE__ << " : " << std::endl 01915 << " For vertex at (" << (a->first >> 2) << ", " << (a->first & 3) << ") :" << std::endl 01916 << " Processor " << proc1 << " has " << all_data[2*idx1+1] << " as the owning part" << std::endl 01917 << " Processor " << proc2 << " has " << all_data[2*idx2+1] << " as the owning part" << std::endl; 01918 } 01919 } 01920 } 01921 return ierr; 01922 } 01923 01924 static int get_part_boundary_verts( iMesh_Instance imesh, 01925 iMeshP_PartitionHandle prtn, 01926 const PartMap& map, 01927 iMeshP_PartHandle part, 01928 std::vector<iBase_EntityHandle>& boundary ) 01929 { 01930 int ierr, logical_id; 01931 ierr = map.part_from_coords( imesh, part, logical_id ); 01932 CHKERR; 01933 01934 int neighbors[5], num_neighbors; 01935 get_part_neighbors( logical_id, map.get_parts().size(), neighbors, num_neighbors ); 01936 01937 for (int j = 0; j < num_neighbors; ++j) { 01938 std::vector<iBase_EntityHandle> iface; 01939 ierr = interface_verts( imesh, prtn, part, neighbors[j], map, iface ); 01940 CHKERR; 01941 std::copy( iface.begin(), iface.end(), std::back_inserter(boundary) ); 01942 } 01943 01944 std::sort( boundary.begin(), boundary.end() ); 01945 boundary.erase( std::unique( boundary.begin(), boundary.end() ), boundary.end() ); 01946 return iBase_SUCCESS; 01947 } 01948 01955 int test_entity_status( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 01956 { 01957 int ierr, rank, size; 01958 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01959 MPI_Comm_size( MPI_COMM_WORLD, &size ); 01960 01961 // get local part handles 01962 std::vector<iMeshP_PartHandle> parts; 01963 ierr = get_local_parts( imesh, prtn, parts ); 01964 PCHECK; 01965 01966 // for each part 01967 int num_quad_ent_incorrect = 0, num_quad_ent_error = 0; 01968 int num_quad_arr_incorrect = 0, num_quad_arr_error = 0; 01969 int num_vert_ent_incorrect = 0, num_vert_ent_error = 0; 01970 int num_vert_arr_incorrect = 0, num_vert_arr_error = 0; 01971 for (size_t i = 0; i < parts.size(); ++i) { 01972 const iMeshP_PartHandle part = parts[i]; 01973 01974 // get quads and vertices 01975 std::vector<iBase_EntityHandle> quads, verts; 01976 ierr = get_part_quads_and_verts( imesh, part, quads, verts ); 01977 if (ierr) 01978 break; 01979 01980 // check quad status (no ghosting yet) 01981 for (size_t j = 0; j < quads.size(); ++j) { 01982 int status; 01983 iMeshP_getEntStatus( imesh, prtn, part, quads[j], &status, &ierr ); 01984 if (ierr != iBase_SUCCESS) { 01985 ++num_quad_ent_error; 01986 ierr = iBase_SUCCESS; 01987 continue; 01988 } 01989 01990 if (status != iMeshP_INTERNAL) 01991 ++num_quad_ent_incorrect; 01992 } 01993 01994 // check quad status using iMeshP_getEntStatusArr 01995 std::vector<int> stat_list(quads.size()); 01996 int* junk1 = &stat_list[0]; 01997 int junk2 = stat_list.size(), count; 01998 iMeshP_getEntStatusArr( imesh, prtn, part, &quads[0], quads.size(), 01999 &junk1, &junk2, &count, &ierr ); 02000 if (ierr != iBase_SUCCESS) { 02001 ++num_quad_arr_error; 02002 ierr = iBase_SUCCESS; 02003 continue; 02004 } 02005 assert( junk1 == &stat_list[0] ); 02006 assert( junk2 == (int)stat_list.size() ); 02007 assert( count == (int)quads.size() ); 02008 for (size_t j = 0; j < quads.size(); ++j) 02009 if (stat_list[j] != iMeshP_INTERNAL) 02010 ++num_quad_arr_incorrect; 02011 02012 // figure out which vertices are on the boundary 02013 std::vector<iBase_EntityHandle> boundary; 02014 ierr = get_part_boundary_verts(imesh, prtn, map, part, boundary); 02015 if (ierr) 02016 break; 02017 std::sort( boundary.begin(), boundary.end() ); 02018 02019 // check vertex status (no ghosting yet) 02020 for (size_t j = 0; j < verts.size(); ++j) { 02021 int status; 02022 iMeshP_getEntStatus( imesh, prtn, part, verts[j], &status, &ierr ); 02023 if (ierr != iBase_SUCCESS) { 02024 ++num_vert_ent_error; 02025 ierr = iBase_SUCCESS; 02026 continue; 02027 } 02028 bool on_boundary = std::binary_search( boundary.begin(), boundary.end(), verts[j] ); 02029 if (status != (on_boundary ? iMeshP_BOUNDARY : iMeshP_INTERNAL)) 02030 ++num_vert_ent_incorrect; 02031 } 02032 02033 // check vert status using iMeshP_getEntStatusArr 02034 stat_list.resize(verts.size()); 02035 junk1 = &stat_list[0]; 02036 junk2 = stat_list.size(); 02037 iMeshP_getEntStatusArr( imesh, prtn, part, &verts[0], verts.size(), 02038 &junk1, &junk2, &count, &ierr ); 02039 if (ierr != iBase_SUCCESS) { 02040 ++num_vert_arr_error; 02041 ierr = iBase_SUCCESS; 02042 continue; 02043 } 02044 assert( junk1 == &stat_list[0] ); 02045 assert( junk2 == (int)stat_list.size() ); 02046 assert( count == (int)verts.size() ); 02047 for (size_t j = 0; j < verts.size(); ++j) { 02048 bool on_boundary = std::binary_search( boundary.begin(), boundary.end(), verts[j] ); 02049 if (stat_list[j] != (on_boundary ? iMeshP_BOUNDARY : iMeshP_INTERNAL)) 02050 ++num_vert_arr_incorrect; 02051 } 02052 } 02053 PCHECK; // check if loop interrupted by any internal errors 02054 02055 ASSERT( 0 == num_quad_ent_error ); 02056 ASSERT( 0 == num_quad_arr_error ); 02057 ASSERT( 0 == num_vert_ent_error ); 02058 ASSERT( 0 == num_vert_arr_error ); 02059 ASSERT( 0 == num_quad_ent_incorrect ); 02060 ASSERT( 0 == num_quad_arr_incorrect ); 02061 ASSERT( 0 == num_vert_ent_incorrect ); 02062 ASSERT( 0 == num_vert_arr_incorrect ); 02063 02064 return iBase_SUCCESS; 02065 } 02066 02073 int test_entity_copy_parts( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 02074 { 02075 int ierr, rank, size; 02076 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02077 MPI_Comm_size( MPI_COMM_WORLD, &size ); 02078 02079 // get local part handles 02080 std::vector<iMeshP_PartHandle> parts; 02081 ierr = get_local_parts( imesh, prtn, parts ); 02082 PCHECK; 02083 ASSERT( !parts.empty() ); 02084 02085 // select a singe part to test 02086 const iMeshP_PartHandle part = parts[0]; 02087 int logical_id; 02088 ierr = map.part_from_coords( imesh, part, logical_id ); 02089 CHKERR; 02090 const iMeshP_Part part_id = map.part_id_from_local_id( logical_id ); 02091 02092 // get vertices in part 02093 std::vector<iBase_EntityHandle> quads, verts; 02094 ierr = get_part_quads_and_verts( imesh, part, quads, verts ); 02095 PCHECK; 02096 02097 // get neighbors 02098 int neighbors[5], num_neighbors; 02099 get_part_neighbors( logical_id, map.get_parts().size(), neighbors, num_neighbors ); 02100 02101 // build map of sharing data for each vertex 02102 std::map< iBase_EntityHandle, std::vector<iMeshP_Part> > vert_sharing; 02103 for (int j = 0; j < num_neighbors; ++j) { 02104 std::vector<iBase_EntityHandle> iface; 02105 ierr = interface_verts( imesh, prtn, part, neighbors[j], map, iface ); 02106 CHKERR; 02107 for (size_t k = 0; k < iface.size(); ++k) 02108 vert_sharing[iface[k]].push_back( map.part_id_from_local_id( neighbors[j] ) ); 02109 } 02110 02111 // test getNumCopies for each vertex 02112 std::map< iBase_EntityHandle, std::vector<iMeshP_Part> >::iterator i; 02113 int num_failed = 0, num_incorrect = 0; 02114 for (i = vert_sharing.begin(); i != vert_sharing.end(); ++i) { 02115 int count; 02116 iBase_EntityHandle vtx = i->first; 02117 iMeshP_getNumCopies( imesh, prtn, vtx, &count, &ierr ); 02118 if (ierr) 02119 ++num_failed; 02120 else if ((unsigned)count != i->second.size()+1) // add one for the part we queried from 02121 ++num_incorrect; 02122 } 02123 ASSERT( 0 == num_failed ); 02124 ASSERT( 0 == num_incorrect ); 02125 02126 // get getCopyParts for each vertex 02127 num_failed = num_incorrect = 0; 02128 for (i = vert_sharing.begin(); i != vert_sharing.end(); ++i) { 02129 iMeshP_Part* list = 0; 02130 int junk = 0, count; 02131 iMeshP_getCopyParts( imesh, prtn, i->first, &list, &junk, &count, &ierr ); 02132 if (iBase_SUCCESS != ierr) { 02133 ++num_failed; 02134 continue; 02135 } 02136 if ((unsigned)count != i->second.size()+1) { // add one for the part we queried from 02137 ++num_incorrect; 02138 free(list); 02139 continue; 02140 } 02141 02142 std::vector<iMeshP_Part> expected( i->second ); 02143 expected.push_back( part_id ); 02144 std::sort( list, list+count ); 02145 std::sort( expected.begin(), expected.end() ); 02146 bool eq = std::equal( list, list+count, expected.begin() ); 02147 free( list ); 02148 if (!eq) 02149 ++num_incorrect; 02150 } 02151 ASSERT( 0 == num_failed ); 02152 ASSERT( 0 == num_incorrect ); 02153 02154 return iBase_SUCCESS; 02155 } 02156 02157 // store remote handle data for a vertex 02158 struct VtxCopyData { 02159 std::vector<iMeshP_Part> parts; 02160 std::vector<iBase_EntityHandle> handles; 02161 }; 02162 02170 int test_entity_copies( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& /* map */ ) 02171 { 02172 int ierr, rank, size; 02173 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02174 MPI_Comm_size( MPI_COMM_WORLD, &size ); 02175 02176 // generate a unique ID for each vertex using the coordinates. 02177 // see create_mesh(..): each vertex has integer coordinates (x,y,0) 02178 // with x in [0,inf] and y in [0,4] 02179 // then to an Allgatherv to exchange handles for each processor 02180 02181 // cast everything to iBase_EntityHandle so we can pack it all in one communication 02182 MPI_Datatype tmp_type; 02183 if (sizeof(iBase_EntityHandle) == sizeof(unsigned)) 02184 tmp_type = MPI_UNSIGNED; 02185 else if (sizeof(iBase_EntityHandle) == sizeof(unsigned long)) 02186 tmp_type = MPI_UNSIGNED_LONG; 02187 else 02188 return iBase_FAILURE; 02189 const MPI_Datatype type = tmp_type; // make it const 02190 02191 // get local part handles 02192 std::vector<iMeshP_PartHandle> parts; 02193 ierr = get_local_parts( imesh, prtn, parts ); 02194 PCHECK; 02195 std::vector<iMeshP_Part> part_ids(parts.size()); 02196 iMeshP_Part* junk1 = &part_ids[0]; 02197 int junk2 = part_ids.size(), junk3; 02198 iMeshP_getPartIdsFromPartHandlesArr( imesh, prtn, &parts[0], parts.size(), 02199 &junk1, &junk2, &junk3, &ierr ); 02200 PCHECK; 02201 assert( junk1 == &part_ids[0] ); 02202 assert( junk2 == (int)part_ids.size() ); 02203 assert( junk3 == (int)parts.size() ); 02204 02205 // build list of {vtx_id, part_id, handle} tuples to send 02206 // also build list of local vertex handles 02207 std::vector<iBase_EntityHandle> local_data, local_vertices; 02208 for (size_t i = 0; i < parts.size(); ++i) { 02209 // get vertices 02210 std::vector<iBase_EntityHandle> quads, verts; 02211 ierr = get_part_quads_and_verts( imesh, parts[i], quads, verts ); 02212 if (ierr) 02213 break; 02214 02215 // add all vertices to local_data 02216 for (size_t j = 0; j < verts.size(); ++j) { 02217 int tag; 02218 ierr = vertex_tag( imesh, verts[j], tag ); 02219 if (ierr) 02220 break; 02221 long tmp_h = tag; 02222 local_data.push_back((iBase_EntityHandle)tmp_h); 02223 tmp_h = part_ids[i]; 02224 local_data.push_back((iBase_EntityHandle)tmp_h); 02225 local_data.push_back( verts[j] ); 02226 } 02227 if (ierr) 02228 break; 02229 02230 std::copy( verts.begin(), verts.end(), std::back_inserter(local_vertices) ); 02231 } 02232 02233 // build list of local vertices 02234 std::sort( local_vertices.begin(), local_vertices.end() ); 02235 local_vertices.erase( std::unique( local_vertices.begin(), local_vertices.end() ), local_vertices.end() ); 02236 std::vector<int> local_vtx_tags(local_vertices.size()); 02237 CHKERR; 02238 for (size_t i = 0; i < local_vertices.size(); ++i) { 02239 ierr = vertex_tag( imesh, local_vertices[i], local_vtx_tags[i] ); 02240 if (ierr) 02241 break; 02242 } 02243 CHKERR; 02244 02245 // communicate data 02246 std::vector<int> gcounts(size), gdisp(size); 02247 int local_data_size = local_data.size(); 02248 ierr = MPI_Allgather( &local_data_size, 1, MPI_INT, &gcounts[0], 1, MPI_INT, MPI_COMM_WORLD ); 02249 CHKERR; 02250 gdisp[0] = 0; 02251 for (int i = 1; i < size; ++i) 02252 gdisp[i] = gdisp[i-1]+gcounts[i-1]; 02253 std::vector<iBase_EntityHandle> global_data( gdisp[size-1]+gcounts[size-1] ); 02254 ierr = MPI_Allgatherv( &local_data[0], local_data_size, type, 02255 &global_data[0], &gcounts[0], &gdisp[0], type, MPI_COMM_WORLD ); 02256 CHKERR; 02257 02258 // arrange global data in a more useful way 02259 std::map<int,VtxCopyData> vtx_sharing; 02260 assert( global_data.size() % 3 == 0 ); 02261 for (size_t i = 0; i < global_data.size(); i += 3) { 02262 int tag = (int)(size_t)global_data[i ]; 02263 iMeshP_Part part = (iMeshP_Part)(size_t)global_data[i+1]; 02264 iBase_EntityHandle handle = global_data[i+2]; 02265 vtx_sharing[tag].parts.push_back( part ); 02266 vtx_sharing[tag].handles.push_back( handle ); 02267 } 02268 02269 // test iMeshP_getCopies for each local vertex 02270 int num_error = 0, num_incorrect = 0, junk4; 02271 for (size_t i = 0; i < local_vertices.size(); ++i) { 02272 int num_copies = -1; 02273 //iMeshP_Part* part_ids = 0; 02274 iMeshP_Part* ptr_part_ids = 0; // Use ptr_part_ids to avoid shadowing std::vector<iMeshP_Part> part_ids 02275 iBase_EntityHandle* copies = 0; 02276 junk2 = junk3 = junk4 = 0; 02277 iMeshP_getCopies( imesh, prtn, local_vertices[i], 02278 &ptr_part_ids, &junk2, &num_copies, 02279 &copies, &junk3, &junk4, &ierr ); 02280 if (iBase_SUCCESS != ierr) { 02281 ++num_error; 02282 continue; 02283 } 02284 assert( junk4 == num_copies ); 02285 02286 VtxCopyData& expected = vtx_sharing[local_vtx_tags[i]]; 02287 if (num_copies != (int)expected.parts.size()) 02288 ++num_incorrect; 02289 else for (size_t j = 0; j < expected.parts.size(); ++j) { 02290 int idx = std::find( ptr_part_ids, ptr_part_ids + num_copies, expected.parts[j] ) - ptr_part_ids; 02291 if (idx == num_copies || copies[idx] != expected.handles[j]) { 02292 ++num_incorrect; 02293 break; 02294 } 02295 } 02296 free(ptr_part_ids); 02297 free(copies); 02298 } 02299 ASSERT( 0 == num_error ); 02300 ASSERT( 0 == num_incorrect ); 02301 02302 // test iMeshP_getCopyOnPart for each local vertex 02303 num_error = num_incorrect = 0; 02304 for (size_t i = 0; i < local_vertices.size(); ++i) { 02305 VtxCopyData& expected = vtx_sharing[local_vtx_tags[i]]; 02306 for (size_t j = 0; j < expected.parts.size(); ++j) { 02307 iBase_EntityHandle copy; 02308 iMeshP_getCopyOnPart( imesh, prtn, local_vertices[i], expected.parts[j], ©, &ierr ); 02309 if (iBase_SUCCESS != ierr) 02310 ++num_error; 02311 else if (expected.handles[j] != copy) 02312 ++num_incorrect; 02313 } 02314 } 02315 ASSERT( 0 == num_error ); 02316 ASSERT( 0 == num_incorrect ); 02317 02318 // test iMeshP_getOwnerCopy for each local vertex 02319 num_error = num_incorrect = 0; 02320 for (size_t i = 0; i < local_vertices.size(); ++i) { 02321 VtxCopyData& expected = vtx_sharing[local_vtx_tags[i]]; 02322 iMeshP_Part owner_id = 0; 02323 iMeshP_getEntOwnerPart( imesh, prtn, local_vertices[i], &owner_id, &ierr ); 02324 if (iBase_SUCCESS != ierr) 02325 continue; // not testing getEntOwnerPart here 02326 02327 size_t idx = std::find( expected.parts.begin(), expected.parts.end(), owner_id ) 02328 - expected.parts.begin(); 02329 if (idx == expected.parts.size()) 02330 continue; // not testing getEntOwnerPart here 02331 02332 iMeshP_Part owner_id_2 = 0; 02333 iBase_EntityHandle copy = 0; 02334 iMeshP_getOwnerCopy( imesh, prtn, local_vertices[i], &owner_id_2, ©, &ierr ); 02335 if (iBase_SUCCESS != ierr) 02336 ++num_error; 02337 else if (owner_id_2 != owner_id && copy != expected.handles[idx]) 02338 ++num_incorrect; 02339 } 02340 ASSERT( 0 == num_error ); 02341 ASSERT( 0 == num_incorrect ); 02342 02343 return iBase_SUCCESS; 02344 } 02345 02346 02347 int get_num_adj_quads( iMesh_Instance imesh, iBase_EntityHandle vtx, int& num ) 02348 { 02349 iBase_EntityHandle* list = 0; 02350 int ierr, junk = 0; 02351 iMesh_getEntAdj( imesh, vtx, iBase_FACE, &list, &junk, &num, &ierr ); 02352 if (iBase_SUCCESS == ierr) 02353 free(list); 02354 return ierr; 02355 } 02356 02357 int get_adj( iMesh_Instance imesh, 02358 iBase_EntityHandle ent, 02359 int type, 02360 std::vector<iBase_EntityHandle>& adj ) 02361 { 02362 iBase_EntityHandle* list = 0; 02363 int ierr, num, junk = 0; 02364 iMesh_getEntAdj( imesh, ent, type, &list, &junk, &num, &ierr ); 02365 if (iBase_SUCCESS == ierr) { 02366 std::copy( list, list+num, std::back_inserter(adj) ); 02367 free( list ); 02368 } 02369 return ierr; 02370 } 02371 02372 // assume regular quad mesh 02373 int get_boundary_vertices( iMesh_Instance imesh, std::vector<iBase_EntityHandle>& bdry ) 02374 { 02375 int ierr, n; 02376 iBase_EntitySetHandle root; 02377 iMesh_getRootSet( imesh, &root, &ierr ); 02378 CHKERR; 02379 std::vector<iBase_EntityHandle> all_verts; 02380 ierr = get_entities( imesh, root, iBase_VERTEX, iMesh_POINT, all_verts ); 02381 CHKERR; 02382 bdry.clear(); 02383 for (size_t i = 0; i < all_verts.size(); ++i) { 02384 ierr = get_num_adj_quads( imesh, all_verts[i], n ); 02385 CHKERR; 02386 if (n != 4) 02387 bdry.push_back( all_verts[i] ); 02388 } 02389 return iBase_SUCCESS; 02390 } 02391 02392 int check_one_layer( iMesh_Instance imesh, iBase_EntityHandle vtx, 02393 const std::vector<iBase_EntityHandle>& sorted_vertices ) 02394 { 02395 int ierr; 02396 if (std::binary_search( sorted_vertices.begin(), sorted_vertices.end(), vtx )) 02397 return iBase_SUCCESS; 02398 std::vector<iBase_EntityHandle> quads, verts; 02399 ierr = get_adj( imesh, vtx, iBase_FACE, quads ); 02400 CHKERR; 02401 for (size_t i = 0; i < quads.size(); ++i) { 02402 verts.clear(); 02403 ierr = get_adj( imesh, quads[i], iBase_VERTEX, verts ); 02404 CHKERR; 02405 for (size_t j = 0; j < verts.size(); ++j) { 02406 if (std::binary_search( sorted_vertices.begin(), sorted_vertices.end(), verts[j] )) 02407 return iBase_SUCCESS; 02408 } 02409 } 02410 02411 return iBase_FAILURE; 02412 } 02413 02414 // get number of adjacent quads to each vertex, both on the local 02415 // processor and in the entire mesh 02416 int get_num_adj_all( iMesh_Instance imesh, 02417 const std::vector<iBase_EntityHandle>& verts, 02418 std::vector<int>& num_local_adj, 02419 std::vector<int>& num_all_adj ) 02420 { 02421 int ierr, size; 02422 MPI_Comm_size( MPI_COMM_WORLD, &size ); 02423 02424 std::vector<int> vtx_tags(verts.size()); 02425 num_local_adj.resize( verts.size() ); 02426 for (size_t i = 0; i < verts.size(); ++i) { 02427 ierr = get_num_adj_quads( imesh, verts[i], num_local_adj[i] ); 02428 CHKERR; 02429 ierr = vertex_tag( imesh, verts[i], vtx_tags[i] ); 02430 CHKERR; 02431 } 02432 02433 std::vector<int> counts(size), displ(size); 02434 int num_vtx = verts.size(); 02435 ierr = MPI_Allgather( &num_vtx, 1, MPI_INT, &counts[0], 1, MPI_INT, MPI_COMM_WORLD ); 02436 CHKERR; 02437 displ[0] = 0; 02438 for (int i = 1; i < size; ++i) 02439 displ[i] = displ[i-1]+counts[i-1]; 02440 int total = displ[size-1]+counts[size-1]; 02441 std::vector<int> all_tags(total), all_adj_counts(total); 02442 ierr = MPI_Allgatherv( &vtx_tags[0], vtx_tags.size(), MPI_INT, &all_tags[0], &counts[0], &displ[0], MPI_INT, MPI_COMM_WORLD ); 02443 CHKERR; 02444 ierr = MPI_Allgatherv( &num_local_adj[0], num_local_adj.size(), MPI_INT, &all_adj_counts[0], &counts[0], &displ[0], MPI_INT, MPI_COMM_WORLD ); 02445 CHKERR; 02446 02447 num_all_adj.clear(); 02448 num_all_adj.resize(total,0); 02449 for (int i = 0; i < total; ++i) { 02450 std::vector<int>::iterator it = std::find( vtx_tags.begin(), vtx_tags.end(), all_tags[i] ); 02451 if (it == vtx_tags.end()) 02452 continue; 02453 int idx = it - vtx_tags.begin(); 02454 num_all_adj[idx] += all_adj_counts[i]; 02455 } 02456 02457 return iBase_SUCCESS; 02458 } 02459 02460 02466 int test_create_ghost_ents( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& /* map */) 02467 { 02468 int ierr; 02469 02470 // get boundary vertices 02471 std::vector<iBase_EntityHandle> bdry; 02472 ierr = get_boundary_vertices( imesh, bdry ); 02473 PCHECK; 02474 // get counts of adjacent entities 02475 std::vector<int> num_local_adj, num_global_adj; 02476 ierr = get_num_adj_all( imesh, bdry, num_local_adj, num_global_adj ); 02477 PCHECK; 02478 // create one layer of ghost entities 02479 iMeshP_createGhostEntsAll( imesh, prtn, iBase_FACE, iBase_VERTEX, 1, 0, &ierr ); 02480 PCHECK; 02481 // check that each vertex has the correct number of adjacent entities 02482 int num_incorrect = 0; 02483 for (size_t i = 0; i < bdry.size(); ++i) { 02484 int n; 02485 ierr = get_num_adj_quads( imesh, bdry[i], n ); 02486 if (iBase_SUCCESS != ierr || num_global_adj[i] != n) 02487 ++num_incorrect; 02488 } 02489 ASSERT( 0 == num_incorrect ); 02490 // get new the new boundary 02491 std::vector<iBase_EntityHandle> new_bdry; 02492 ierr = get_boundary_vertices( imesh, new_bdry ); 02493 PCHECK; 02494 // check that each vertex on the new boundary is separated by 02495 // at most one layer from the old boundary 02496 std::sort( bdry.begin(), bdry.end() ); 02497 num_incorrect = 0; 02498 for (size_t i = 0; i < new_bdry.size(); ++i) { 02499 ierr = check_one_layer( imesh, new_bdry[i], bdry ); 02500 if (ierr) 02501 ++num_incorrect; 02502 } 02503 ASSERT( 0 == num_incorrect ); 02504 // make another layer of ghost entiites 02505 bdry.swap( new_bdry ); 02506 new_bdry.clear(); 02507 ierr = get_num_adj_all( imesh, bdry, num_local_adj, num_global_adj ); 02508 PCHECK; 02509 iMeshP_createGhostEntsAll( imesh, prtn, iBase_FACE, iBase_VERTEX, 2, 0, &ierr ); 02510 PCHECK; 02511 // check that each vertex has the correct number of adjacent entities 02512 num_incorrect = 0; 02513 for (size_t i = 0; i < bdry.size(); ++i) { 02514 int n; 02515 ierr = get_num_adj_quads( imesh, bdry[i], n ); 02516 if (iBase_SUCCESS != ierr || num_global_adj[i] != n) 02517 ++num_incorrect; 02518 } 02519 // check that each vertex on the new boundary is separated by 02520 // at most one layer from the old boundary 02521 std::sort( bdry.begin(), bdry.end() ); 02522 num_incorrect = 0; 02523 for (size_t i = 0; i < new_bdry.size(); ++i) { 02524 ierr = check_one_layer( imesh, new_bdry[i], bdry ); 02525 if (ierr) 02526 ++num_incorrect; 02527 } 02528 ASSERT( 0 == num_incorrect ); 02529 02530 return iBase_SUCCESS; 02531 } 02532 02538 int test_exchange_ents( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& map ) 02539 { 02540 int ierr, rank, size; 02541 int num_err = 0; 02542 iMeshP_RequestHandle request; 02543 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02544 MPI_Comm_size( MPI_COMM_WORLD, &size ); 02545 02546 std::vector<iBase_EntityHandle> all_elems; 02547 std::vector<iMeshP_Part> all_ids; 02548 std::vector<iBase_EntityHandle> quads; 02549 02550 // get local part handles and part ids 02551 std::vector<iMeshP_PartHandle> local_handles; 02552 std::vector<iMeshP_Part> local_ids; 02553 ierr = get_local_parts( imesh, prtn, local_handles, &local_ids ); 02554 PCHECK; 02555 02556 // get loacal quads before exchange 02557 quads.clear(); 02558 ierr = get_entities( imesh, local_handles[0], iBase_FACE, iMesh_QUADRILATERAL, quads ); 02559 CHKERR; 02560 int n_quads = quads.size(); 02561 02562 // send all elements in local processor to all other processors 02563 for (size_t i = 0; i < map.get_parts().size(); ++i) { 02564 if (map.get_parts()[i] == (unsigned int) rank) continue; // skip own rank 02565 02566 for (int j = 0; j < n_quads; j++) { 02567 all_elems.push_back(quads[j]); 02568 all_ids.push_back(map.get_parts()[i]); 02569 } 02570 } 02571 02572 // exchange entities 02573 iMeshP_exchEntArrToPartsAll(imesh, prtn, &all_elems[0], all_elems.size(), 02574 &all_ids[0], 0, 0, &request, &ierr); 02575 if (iBase_SUCCESS != ierr) ++num_err; 02576 02577 // get local quads after exchange 02578 quads.clear(); 02579 ierr = get_entities( imesh, local_handles[0], iBase_FACE, iMesh_QUADRILATERAL, quads ); 02580 CHKERR; 02581 02582 // # of elements should be # of quads * # of processors 02583 ASSERT(quads.size() == (unsigned int) n_quads*size); 02584 02585 ASSERT(0 == num_err); 02586 02587 return iBase_SUCCESS; 02588 } 02589 02596 int test_push_tag_data_common( iMesh_Instance imesh, 02597 iMeshP_PartitionHandle prtn, 02598 int num_ghost_layers ) 02599 { 02600 const char* src_name = "test_src"; 02601 const char* dst_name = "test_dst"; 02602 int ierr, rank; 02603 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02604 02605 if (num_ghost_layers) { 02606 iMeshP_createGhostEntsAll( imesh, prtn, iBase_FACE, iBase_VERTEX, num_ghost_layers, 0, &ierr ); 02607 PCHECK; 02608 } 02609 02610 iBase_TagHandle src_tag, dst_tag; 02611 iMesh_createTag( imesh, src_name, 1, iBase_INTEGER, &src_tag, &ierr, strlen(src_name) ); 02612 CHKERR; 02613 iMesh_createTag( imesh, dst_name, 1, iBase_INTEGER, &dst_tag, &ierr, strlen(dst_name) ); 02614 CHKERR; 02615 02616 iBase_EntitySetHandle root; 02617 iMesh_getRootSet( imesh, &root, &ierr ); CHKERR; 02618 02619 std::vector<iBase_EntityHandle> verts; 02620 ierr = get_entities( imesh, root, iBase_VERTEX, iMesh_POINT, verts ); 02621 CHKERR; 02622 02623 // test iMeshP_pushTags 02624 // each processor writes its rank on all vertices 02625 // after push, each vertex should be tagged with the rank of its owner 02626 02627 std::vector<int> tag_vals( verts.size(), rank ); 02628 iMesh_setIntArrData( imesh, &verts[0], verts.size(), src_tag, &tag_vals[0], tag_vals.size(), &ierr ); 02629 CHKERR; 02630 02631 iMeshP_pushTags( imesh, prtn, src_tag, dst_tag, iBase_VERTEX, iMesh_POINT, &ierr ); 02632 PCHECK; 02633 02634 tag_vals.clear(); 02635 tag_vals.resize( verts.size(), -1 ); 02636 iBase_TagHandle id_tag; 02637 iMesh_getTagHandle( imesh, "GLOBAL_ID", &id_tag, &ierr, strlen("GLOBAL_ID") ); 02638 std::vector<int> ids(verts.size()); 02639 int* junk1 = &ids[0], junk2 = ids.size(), junk3; 02640 iMesh_getIntArrData( imesh, &verts[0], verts.size(), id_tag, &junk1, &junk2, &junk3, &ierr ); 02641 PCHECK; 02642 int errcount = 0; 02643 for (size_t i = 0; i < verts.size(); ++i) { 02644 iMesh_getIntData( imesh, verts[i], dst_tag, &tag_vals[i], &ierr ); 02645 if (ierr != iBase_SUCCESS) { 02646 std::cerr << "Rank " << rank << " : getIntData failed for vertex " << ids[i] << std::endl; 02647 std::cerr.flush(); 02648 ++errcount; 02649 } 02650 } 02651 ASSERT(0 == errcount); 02652 02653 // int *junk1 = &tag_vals[0], junk2 = tag_vals.size(), junk3; 02654 // iMesh_getIntArrData( imesh, &verts[0], verts.size(), dst_tag, &junk1, &junk2, &junk3, &ierr ); 02655 // PCHECK; 02656 // assert( junk1 == &tag_vals[0] ); 02657 // assert( junk2 == (int)tag_vals.size() ); 02658 // assert( junk3 == (int)verts.size() ); 02659 02660 std::vector<int> expected( verts.size() ); 02661 std::vector<iMeshP_Part> parts( verts.size() ); 02662 iMeshP_Part* junk4 = &parts[0]; 02663 junk2 = parts.size(); 02664 iMeshP_getEntOwnerPartArr( imesh, prtn, &verts[0], verts.size(), &junk4, &junk2, &junk3, &ierr ); 02665 PCHECK; 02666 assert(junk4 == &parts[0]); 02667 assert(junk2 == (int)parts.size()); 02668 assert(junk3 == (int)verts.size()); 02669 junk1 = &expected[0]; 02670 junk2 = expected.size(); 02671 iMeshP_getRankOfPartArr( imesh, prtn, &parts[0], parts.size(), &junk1, &junk2, &junk3, &ierr ); 02672 PCHECK; 02673 assert(junk1 == &expected[0]); 02674 assert(junk2 == (int)expected.size()); 02675 assert(junk3 == (int)parts.size()); 02676 02677 ASSERT( tag_vals == expected ); 02678 02679 02680 02681 // test iMeshP_pushTagsEnt 02682 // write -1 on all vertices 02683 // For each vertex owned by this processor and shared with more than 02684 // two others, write the rank of the owning processor. 02685 02686 tag_vals.clear(); 02687 tag_vals.resize( verts.size(), -1 ); 02688 iMesh_setIntArrData( imesh, &verts[0], verts.size(), src_tag, &tag_vals[0], tag_vals.size(), &ierr ); 02689 PCHECK; 02690 tag_vals.resize( verts.size(), -1 ); 02691 iMesh_setIntArrData( imesh, &verts[0], verts.size(), dst_tag, &tag_vals[0], tag_vals.size(), &ierr ); 02692 PCHECK; 02693 02694 std::vector<iBase_EntityHandle> some; 02695 for (size_t i = 0; i < verts.size(); ++i) { 02696 int num; 02697 iMeshP_getNumCopies( imesh, prtn, verts[i], &num, &ierr ); 02698 if (iBase_SUCCESS != ierr) 02699 break; 02700 if (num > 2) 02701 some.push_back( verts[i] ); 02702 else 02703 expected[i] = -1; 02704 } 02705 02706 tag_vals.clear(); 02707 tag_vals.resize( some.size(), rank ); 02708 iMesh_setIntArrData( imesh, &some[0], some.size(), src_tag, &tag_vals[0], tag_vals.size(), &ierr ); 02709 PCHECK; 02710 02711 iMeshP_pushTagsEnt( imesh, prtn, src_tag, dst_tag, &some[0], some.size(), &ierr ); 02712 PCHECK; 02713 02714 tag_vals.clear(); 02715 tag_vals.resize( verts.size(), -1 ); 02716 junk1 = &tag_vals[0]; 02717 junk2 = tag_vals.size(); 02718 iMesh_getIntArrData( imesh, &verts[0], verts.size(), dst_tag, &junk1, &junk2, &junk3, &ierr ); 02719 CHKERR; 02720 assert( junk1 == &tag_vals[0] ); 02721 assert( junk2 == (int)tag_vals.size() ); 02722 assert( junk3 == (int)verts.size() ); 02723 02724 ASSERT( tag_vals == expected ); 02725 return iBase_SUCCESS; 02726 } 02727 02734 int test_push_tag_data_iface( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& ) 02735 { 02736 return test_push_tag_data_common( imesh, prtn, 0 ); 02737 } 02738 02745 int test_push_tag_data_ghost( iMesh_Instance imesh, iMeshP_PartitionHandle prtn, const PartMap& ) 02746 { 02747 return test_push_tag_data_common( imesh, prtn, 1 ); 02748 } 02749 02750 02751 /************************************************************************** 02752 PartMap class 02753 **************************************************************************/ 02754 02755 02756 02757 02758 int PartMap::build_map( iMesh_Instance imesh, 02759 iMeshP_PartitionHandle prtn, 02760 int num_expected_parts ) 02761 { 02762 int ierr, rank, size; 02763 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02764 MPI_Comm_size( MPI_COMM_WORLD, &size ); 02765 02766 // get local parts 02767 std::vector<iMeshP_PartHandle> local_parts; 02768 std::vector<iMeshP_Part> imesh_ids; 02769 ierr = get_local_parts( imesh, prtn, local_parts, &imesh_ids ); 02770 CHKERR; 02771 02772 // get logical ids for local parts 02773 std::vector<int> local_ids(local_parts.size()); 02774 for (size_t i = 0; i < local_parts.size(); ++i) { 02775 ierr = part_from_coords( imesh, local_parts[i], local_ids[i] ); 02776 CHKERR; 02777 } 02778 02779 // get total number of parts 02780 int num_global = 0, num_local = local_parts.size(); 02781 ierr = MPI_Allreduce( &num_local, &num_global, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); 02782 CHKERR; 02783 if (num_global != num_expected_parts) { 02784 std::cerr << "Invalid/unexpected global part count at " __FILE__ ":" 02785 << __LINE__ << " (proc " << rank << "): " << std::endl 02786 << " Expected: " << num_expected_parts << std::endl 02787 << " Actual: " << num_global << std::endl; 02788 return 1; 02789 } 02790 02791 // get counts and displacements for Allgatherv calls 02792 std::vector<int> dspls(size), counts(size); 02793 ierr = MPI_Allgather( &num_local, 1, MPI_INT, &counts[0], 1, MPI_INT, MPI_COMM_WORLD ); 02794 CHKERR; 02795 dspls[0] = 0; 02796 for (int i = 1; i < size; ++i) 02797 dspls[i] = dspls[i-1] + counts[i-1]; 02798 02799 // gather iMeshP_Part list from each processor 02800 std::vector<unsigned> global_part_ids(num_expected_parts); 02801 assert(sizeof(iMeshP_Part) == sizeof(int)); 02802 ierr = MPI_Allgatherv( &imesh_ids[0], num_local, MPI_UNSIGNED, 02803 &global_part_ids[0], &counts[0], &dspls[0], MPI_UNSIGNED, 02804 MPI_COMM_WORLD ); 02805 CHKERR; 02806 02807 // gather local ids from each processor 02808 std::vector<int> global_id_list(num_expected_parts); 02809 ierr = MPI_Allgatherv( &local_ids[0], num_local, MPI_INT, 02810 &global_id_list[0], &counts[0], &dspls[0], MPI_INT, 02811 MPI_COMM_WORLD ); 02812 CHKERR; 02813 02814 // build owner list 02815 std::vector<int> global_owners(num_expected_parts); 02816 for (int i = 0; i < size; ++i) 02817 for (int j = 0; j < counts[i]; ++j) 02818 global_owners[dspls[i]+j] = i; 02819 02820 02821 // populate member lists 02822 sortedPartList = global_part_ids; 02823 std::sort( sortedPartList.begin(), sortedPartList.end() ); 02824 partLocalIds.resize( num_expected_parts ); 02825 partRanks.resize( num_expected_parts ); 02826 for (int i = 0; i < num_expected_parts; ++i) { 02827 int idx = std::lower_bound( sortedPartList.begin(), sortedPartList.end(), global_part_ids[i] ) - sortedPartList.begin(); 02828 partLocalIds[idx] = global_id_list[i]; 02829 partRanks[idx] = global_owners[i]; 02830 } 02831 02832 // do some consistency checking 02833 if (std::unique( sortedPartList.begin(), sortedPartList.end() ) != sortedPartList.end()) { 02834 if (rank == 0) { 02835 std::cerr << "ERROR: Duplicate iMeshP_Part values detected at " __FILE__ ":" << __LINE__ << std::endl; 02836 } 02837 return 1; 02838 } 02839 02840 // build revesre local id map and check for duplicates 02841 localIdReverseMap.clear(); 02842 localIdReverseMap.resize(num_expected_parts, -1); 02843 for (int i = 0; i < num_expected_parts; ++i) { 02844 int idx = partLocalIds[i]; 02845 if (localIdReverseMap[idx] != -1) { 02846 if (rank == 0) { 02847 std::cerr << "ERROR: Part mesh has been duplicated in multiple parts." << std::endl 02848 << " Detected at " __FILE__ ":" << __LINE__ << std::endl 02849 << " See PartMap::part_from_coords" << std::endl; 02850 } 02851 return 1; 02852 } 02853 if (idx >= num_expected_parts) { 02854 if (rank == 0) { 02855 std::cerr << "ERROR: Part mesh invalid/incorrect mesh." << std::endl 02856 << " Detected at " __FILE__ ":" << __LINE__ << std::endl 02857 << " See PartMap::part_from_coords" << std::endl; 02858 } 02859 return 1; 02860 } 02861 02862 localIdReverseMap[idx] = i; 02863 } 02864 02865 return 0; 02866 } 02867 02868 02869 void PartMap::part_id_from_rank( int rank, std::vector<iMeshP_Part>& parts ) const 02870 { 02871 for (size_t i = 0; i < sortedPartList.size(); ++i) 02872 if (partRanks[i] == rank) 02873 parts.push_back( sortedPartList[i] ); 02874 } 02875 02876 void PartMap::local_id_from_rank( int rank, std::vector<int>& ids ) const 02877 { 02878 for (size_t i = 0; i < sortedPartList.size(); ++i) 02879 if (partRanks[i] == rank) 02880 ids.push_back( partLocalIds[i] ); 02881 } 02882 02883 02884 int PartMap::part_from_coords( iMesh_Instance imesh, iMeshP_PartHandle part, int& id ) 02885 { 02886 int ierr, rank; 02887 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 02888 02889 // get elements 02890 const int num_elem = 4; 02891 iBase_EntityHandle array[num_elem]; 02892 iBase_EntityHandle* ptr = array; 02893 int junk1 = num_elem, n = -1; 02894 iMesh_getEntities( imesh, part, iBase_FACE, iMesh_QUADRILATERAL, &ptr, &junk1, &n, &ierr ); 02895 CHKERR; 02896 assert( ptr == array ); 02897 assert( junk1 == num_elem ); 02898 if (n != num_elem) { 02899 std::cerr << "Internal error at " __FILE__ ":" << __LINE__ << " (proc " 02900 << rank << "): Expected all parts to have " << num_elem 02901 << " elements. Found one with " << n << std::endl; 02902 return 1; 02903 } 02904 02905 // get vertices 02906 iBase_EntityHandle adj_array[4*num_elem]; 02907 int junk2, junk3, offset_array[5]; 02908 ptr = adj_array; 02909 junk1 = sizeof(adj_array)/sizeof(adj_array[0]); 02910 junk2 = sizeof(offset_array)/sizeof(offset_array[0]); 02911 int* ptr2 = offset_array; 02912 iMesh_getEntArrAdj( imesh, array, num_elem, iBase_VERTEX, 02913 &ptr, &junk1, &n, 02914 &ptr2, &junk2, &junk3, 02915 &ierr ); 02916 CHKERR; 02917 assert( ptr == adj_array ); 02918 assert( ptr2 == offset_array ); 02919 assert( junk1 == sizeof(adj_array)/sizeof(adj_array[0]) ); 02920 assert( junk2 == sizeof(offset_array)/sizeof(offset_array[0]) ); 02921 assert( n == 4*num_elem ); 02922 assert( offset_array[0] == 0 ); 02923 for (int i = 1; i < junk3; ++i) 02924 assert( offset_array[i] - offset_array[i-1] == 4 ); 02925 02926 // find center vertex 02927 iBase_EntityHandle vtx; 02928 bool all_match; 02929 for (int i = 0; i < 4; ++i) { 02930 vtx = adj_array[i]; 02931 all_match = true; 02932 for (int j = 1; j < 4; ++j) { 02933 iBase_EntityHandle* mvtx = adj_array + 4*j; 02934 int k; 02935 for (k = 0; k < 4; ++k) 02936 if (mvtx[k] == vtx) 02937 break; 02938 if (k == 4) 02939 all_match = false; 02940 } 02941 if (all_match) 02942 break; 02943 } 02944 assert(all_match); 02945 02946 // get center vertex coordinates 02947 double x, y, z; 02948 iMesh_getVtxCoord( imesh, vtx, &x, &y, &z, &ierr ); 02949 CHKERR; 02950 assert( 0.0 == z ); 02951 const int xi = ((int)round(x) - 1)/2; 02952 const int yi = ((int)round(y) - 1)/2; 02953 assert (xi >= 0); 02954 assert (yi >= 0); 02955 assert( fabs(x - 2*xi - 1) < 1e-12 ); 02956 assert( fabs(y - 2*yi - 1) < 1e-12 ); 02957 02958 id = 2*xi + yi; 02959 return 0; 02960 }