MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* $Id$ 00002 * 00003 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00004 * storing and accessing finite element mesh data. 00005 * 00006 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00007 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00008 * retains certain rights in this software. 00009 * 00010 * This library is free software; you can redistribute it and/or 00011 * modify it under the terms of the GNU Lesser General Public 00012 * License as published by the Free Software Foundation; either 00013 * version 2.1 of the License, or (at your option) any later version. 00014 * 00015 */ 00016 00017 /** Get a mesh from MOAB and write a Zoltan partition set back into MOAB and to 00018 * a file 00019 * 00020 * 00021 */ 00022 00023 #include <iostream> 00024 #include <assert.h> 00025 #include <sstream> 00026 #include <algorithm> 00027 00028 #include "moab/ZoltanPartitioner.hpp" 00029 #include "moab/Interface.hpp" 00030 #include "Internals.hpp" 00031 #include "moab/Range.hpp" 00032 #include "moab/WriteUtilIface.hpp" 00033 #include "moab/MeshTopoUtil.hpp" 00034 #include "MBTagConventions.hpp" 00035 #include "moab/CN.hpp" 00036 // used for gnomonic projection 00037 #include "moab/IntxMesh/IntxUtils.hpp" 00038 00039 #ifdef MOAB_HAVE_CGM 00040 #include "CGMConfig.h" 00041 #include <limits> 00042 #include "RefEntity.hpp" 00043 #include "RefFace.hpp" 00044 #include "RefEdge.hpp" 00045 #include "RefVertex.hpp" 00046 #include "Body.hpp" 00047 #include "TopologyEntity.hpp" 00048 #include "CastTo.hpp" 00049 #include "CABodies.hpp" 00050 #include "TDParallel.hpp" 00051 #include "TDUniqueId.hpp" 00052 #endif 00053 00054 using namespace moab; 00055 00056 #define RR \ 00057 if( MB_SUCCESS != result ) return result 00058 00059 static double* Points = NULL; 00060 static int* GlobalIds = NULL; 00061 static int NumPoints = 0; 00062 static int* NumEdges = NULL; 00063 static int* NborGlobalId = NULL; 00064 static int* NborProcs = NULL; 00065 static double* ObjWeights = NULL; 00066 static double* EdgeWeights = NULL; 00067 static int* Parts = NULL; 00068 00069 const bool debug = false; 00070 00071 ZoltanPartitioner::ZoltanPartitioner( Interface* impl, const bool use_coords, int argc, char** argv 00072 #ifdef MOAB_HAVE_CGM 00073 , 00074 GeometryQueryTool* gqt 00075 #endif 00076 ) 00077 : PartitionerBase< int >( impl, use_coords ), myZZ( NULL ), argcArg( argc ), argvArg( argv ) 00078 #ifdef MOAB_HAVE_CGM 00079 , 00080 gti( gqt ) 00081 #endif 00082 { 00083 } 00084 00085 ZoltanPartitioner::~ZoltanPartitioner() 00086 { 00087 if( NULL != myZZ ) delete myZZ; 00088 } 00089 00090 ErrorCode ZoltanPartitioner::balance_mesh( const char* zmethod, const char* other_method, const bool write_as_sets, 00091 const bool write_as_tags ) 00092 { 00093 if( !strcmp( zmethod, "RR" ) && !strcmp( zmethod, "RCB" ) && !strcmp( zmethod, "RIB" ) && 00094 !strcmp( zmethod, "HSFC" ) && !strcmp( zmethod, "Hypergraph" ) && !strcmp( zmethod, "PHG" ) && 00095 !strcmp( zmethod, "PARMETIS" ) && !strcmp( zmethod, "OCTPART" ) ) 00096 { 00097 std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be " 00098 << "RR, RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl; 00099 return MB_FAILURE; 00100 } 00101 00102 std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB 00103 std::vector< int > ids; // point ids from MOAB 00104 std::vector< int > adjs, length; 00105 Range elems; 00106 00107 // Get a mesh from MOAB and divide it across processors. 00108 00109 ErrorCode result; 00110 00111 if( mbpc->proc_config().proc_rank() == 0 ) 00112 { 00113 result = assemble_graph( 3, pts, ids, adjs, length, elems );RR; 00114 } 00115 00116 myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0] ); 00117 00118 // Initialize Zoltan. This is a C call. The simple C++ code 00119 // that creates Zoltan objects does not keep track of whether 00120 // Zoltan_Initialize has been called. 00121 00122 float version; 00123 00124 Zoltan_Initialize( argcArg, argvArg, &version ); 00125 00126 // Create Zoltan object. This calls Zoltan_Create. 00127 if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() ); 00128 00129 if( NULL == zmethod || !strcmp( zmethod, "RCB" ) ) 00130 SetRCB_Parameters(); 00131 else if( !strcmp( zmethod, "RIB" ) ) 00132 SetRIB_Parameters(); 00133 else if( !strcmp( zmethod, "HSFC" ) ) 00134 SetHSFC_Parameters(); 00135 else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) ) 00136 if( NULL == other_method ) 00137 SetHypergraph_Parameters( "auto" ); 00138 else 00139 SetHypergraph_Parameters( other_method ); 00140 else if( !strcmp( zmethod, "PARMETIS" ) ) 00141 { 00142 if( NULL == other_method ) 00143 SetPARMETIS_Parameters( "RepartGDiffusion" ); 00144 else 00145 SetPARMETIS_Parameters( other_method ); 00146 } 00147 else if( !strcmp( zmethod, "OCTPART" ) ) 00148 { 00149 if( NULL == other_method ) 00150 SetOCTPART_Parameters( "2" ); 00151 else 00152 SetOCTPART_Parameters( other_method ); 00153 } 00154 00155 // Call backs: 00156 00157 myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL ); 00158 myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL ); 00159 myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL ); 00160 myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL ); 00161 myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL ); 00162 myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL ); 00163 00164 // Perform the load balancing partitioning 00165 00166 int changes; 00167 int numGidEntries; 00168 int numLidEntries; 00169 int numImport; 00170 ZOLTAN_ID_PTR importGlobalIds; 00171 ZOLTAN_ID_PTR importLocalIds; 00172 int* importProcs; 00173 int* importToPart; 00174 int numExport; 00175 ZOLTAN_ID_PTR exportGlobalIds; 00176 ZOLTAN_ID_PTR exportLocalIds; 00177 int* exportProcs; 00178 int* exportToPart; 00179 00180 int rc = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, numImport, importGlobalIds, importLocalIds, 00181 importProcs, importToPart, numExport, exportGlobalIds, exportLocalIds, exportProcs, 00182 exportToPart ); 00183 00184 rc = mbGlobalSuccess( rc ); 00185 00186 if( !rc ) { mbPrintGlobalResult( "==============Result==============", myNumPts, numImport, numExport, changes ); } 00187 else 00188 { 00189 return MB_FAILURE; 00190 } 00191 00192 // take results & write onto MOAB partition sets 00193 00194 int* assignment; 00195 00196 mbFinalizePoints( (int)ids.size(), numExport, exportLocalIds, exportProcs, &assignment ); 00197 00198 if( mbpc->proc_config().proc_rank() == 0 ) 00199 { 00200 result = write_partition( mbpc->proc_config().proc_size(), elems, assignment, write_as_sets, write_as_tags ); 00201 00202 if( MB_SUCCESS != result ) return result; 00203 00204 free( (int*)assignment ); 00205 } 00206 00207 // Free the memory allocated for lists returned by LB_Parition() 00208 00209 myZZ->LB_Free_Part( &importGlobalIds, &importLocalIds, &importProcs, &importToPart ); 00210 myZZ->LB_Free_Part( &exportGlobalIds, &exportLocalIds, &exportProcs, &exportToPart ); 00211 00212 // Implementation note: A Zoltan object contains an MPI communicator. 00213 // When the Zoltan object is destroyed, it uses it's MPI communicator. 00214 // So it is important that the Zoltan object is destroyed before 00215 // the MPI communicator is destroyed. To ensure this, dynamically 00216 // allocate the Zoltan object, so you can explicitly destroy it. 00217 // If you create a Zoltan object on the stack, it's destructor will 00218 // be invoked atexit, possibly after the communicator's 00219 // destructor. 00220 00221 return MB_SUCCESS; 00222 } 00223 00224 ErrorCode ZoltanPartitioner::repartition( std::vector< double >& x, std::vector< double >& y, std::vector< double >& z, 00225 int StartID, const char* zmethod, Range& localGIDs ) 00226 { 00227 // 00228 int nprocs = mbpc->proc_config().proc_size(); 00229 int rank = mbpc->proc_config().proc_rank(); 00230 clock_t t = clock(); 00231 // form pts and ids, as in assemble_graph 00232 std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB 00233 pts.resize( x.size() * 3 ); 00234 std::vector< int > ids; // point ids from MOAB 00235 ids.resize( x.size() ); 00236 for( size_t i = 0; i < x.size(); i++ ) 00237 { 00238 pts[3 * i] = x[i]; 00239 pts[3 * i + 1] = y[i]; 00240 pts[3 * i + 2] = z[i]; 00241 ids[i] = StartID + (int)i; 00242 } 00243 // do not get out until done! 00244 00245 Points = &pts[0]; 00246 GlobalIds = &ids[0]; 00247 NumPoints = (int)x.size(); 00248 NumEdges = NULL; 00249 NborGlobalId = NULL; 00250 NborProcs = NULL; 00251 ObjWeights = NULL; 00252 EdgeWeights = NULL; 00253 Parts = NULL; 00254 00255 float version; 00256 if( rank == 0 ) std::cout << "Initializing zoltan..." << std::endl; 00257 00258 Zoltan_Initialize( argcArg, argvArg, &version ); 00259 00260 // Create Zoltan object. This calls Zoltan_Create. 00261 if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() ); 00262 00263 if( NULL == zmethod || !strcmp( zmethod, "RCB" ) ) 00264 SetRCB_Parameters(); 00265 else if( !strcmp( zmethod, "RIB" ) ) 00266 SetRIB_Parameters(); 00267 else if( !strcmp( zmethod, "HSFC" ) ) 00268 SetHSFC_Parameters(); 00269 00270 // set # requested partitions 00271 char buff[10]; 00272 sprintf( buff, "%d", nprocs ); 00273 int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff ); 00274 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00275 00276 // request all, import and export 00277 retval = myZZ->Set_Param( "RETURN_LISTS", "ALL" ); 00278 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00279 00280 myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL ); 00281 myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL ); 00282 myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL ); 00283 myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL ); 00284 myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL ); 00285 myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL ); 00286 00287 // Perform the load balancing partitioning 00288 00289 int changes; 00290 int numGidEntries; 00291 int numLidEntries; 00292 int num_import; 00293 ZOLTAN_ID_PTR import_global_ids, import_local_ids; 00294 int* import_procs; 00295 int* import_to_part; 00296 int num_export; 00297 ZOLTAN_ID_PTR export_global_ids, export_local_ids; 00298 int *assign_procs, *assign_parts; 00299 00300 if( rank == 0 ) 00301 std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nprocs 00302 << " processors..." << std::endl; 00303 00304 retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids, 00305 import_procs, import_to_part, num_export, export_global_ids, export_local_ids, 00306 assign_procs, assign_parts ); 00307 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00308 00309 if( rank == 0 ) 00310 { 00311 std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n"; 00312 t = clock(); 00313 } 00314 00315 std::sort( import_global_ids, import_global_ids + num_import, std::greater< int >() ); 00316 std::sort( export_global_ids, export_global_ids + num_export, std::greater< int >() ); 00317 00318 Range iniGids( (EntityHandle)StartID, (EntityHandle)StartID + x.size() - 1 ); 00319 Range imported, exported; 00320 std::copy( import_global_ids, import_global_ids + num_import, range_inserter( imported ) ); 00321 std::copy( export_global_ids, export_global_ids + num_export, range_inserter( exported ) ); 00322 localGIDs = subtract( iniGids, exported ); 00323 localGIDs = unite( localGIDs, imported ); 00324 // Free data structures allocated by Zoltan::LB_Partition 00325 retval = myZZ->LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part ); 00326 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00327 retval = myZZ->LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts ); 00328 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00329 00330 return MB_SUCCESS; 00331 } 00332 00333 ErrorCode ZoltanPartitioner::partition_inferred_mesh( EntityHandle sfileset, size_t num_parts, int part_dim, 00334 const bool write_as_sets, int projection_type ) 00335 { 00336 ErrorCode result; 00337 00338 moab::Range elverts; 00339 result = mbImpl->get_entities_by_dimension( sfileset, part_dim, elverts );RR; 00340 00341 std::vector< double > elcoords( elverts.size() * 3 ); 00342 result = mbImpl->get_coords( elverts, &elcoords[0] );RR; 00343 00344 std::vector< std::vector< EntityHandle > > part_assignments( num_parts ); 00345 int part, proc; 00346 00347 // Loop over coordinates 00348 for( size_t iel = 0; iel < elverts.size(); iel++ ) 00349 { 00350 // Gather coordinates into temporary array 00351 double* ecoords = &elcoords[iel * 3]; 00352 00353 // do a projection if needed 00354 if (projection_type > 0) 00355 IntxUtils::transform_coordinates(ecoords, projection_type); 00356 00357 // Compute the coordinate's part assignment 00358 myZZ->LB_Point_PP_Assign( ecoords, proc, part ); 00359 00360 // Store the part assignment in the return array 00361 part_assignments[part].push_back( elverts[iel] ); 00362 } 00363 00364 // get the partition set tag 00365 Tag part_set_tag = mbpc->partition_tag(); 00366 00367 // If some entity sets exist, let us delete those first. 00368 if( write_as_sets ) 00369 { 00370 moab::Range oldpsets; 00371 result = mbImpl->get_entities_by_type_and_tag( sfileset, MBENTITYSET, &part_set_tag, NULL, 1, oldpsets, 00372 Interface::UNION );RR; 00373 result = mbImpl->remove_entities( sfileset, oldpsets );RR; 00374 } 00375 00376 size_t nparts_assigned = 0; 00377 for( size_t partnum = 0; partnum < num_parts; ++partnum ) 00378 { 00379 std::vector< EntityHandle >& partvec = part_assignments[partnum]; 00380 00381 nparts_assigned += ( partvec.size() ? 1 : 0 ); 00382 00383 if( write_as_sets ) 00384 { 00385 EntityHandle partNset; 00386 result = mbImpl->create_meshset( moab::MESHSET_SET, partNset );RR; 00387 if( partvec.size() ) 00388 { 00389 result = mbImpl->add_entities( partNset, &partvec[0], partvec.size() );RR; 00390 } 00391 result = mbImpl->add_parent_child( sfileset, partNset );RR; 00392 00393 // part numbering should start from 0 00394 int ipartnum = (int)partnum; 00395 00396 // assign partitions as a sparse tag by grouping elements under sets 00397 result = mbImpl->tag_set_data( part_set_tag, &partNset, 1, &ipartnum );RR; 00398 } 00399 else 00400 { 00401 /* assign as a dense vector to all elements 00402 allocate integer-size partitions */ 00403 std::vector< int > assignment( partvec.size(), 00404 partnum ); // initialize all values to partnum 00405 result = mbImpl->tag_set_data( part_set_tag, &partvec[0], partvec.size(), &assignment[0] );RR; 00406 } 00407 } 00408 00409 if( nparts_assigned != num_parts ) 00410 { 00411 std::cout << "WARNING: The inference yielded lesser number of parts (" << nparts_assigned 00412 << ") than requested by user (" << num_parts << ").\n"; 00413 } 00414 00415 return MB_SUCCESS; 00416 } 00417 00418 ErrorCode ZoltanPartitioner::partition_mesh_and_geometry( const double part_geom_mesh_size, const int nparts, 00419 const char* zmethod, const char* other_method, 00420 double imbal_tol, const int part_dim, 00421 const bool write_as_sets, const bool write_as_tags, 00422 const int obj_weight, const int edge_weight, 00423 #ifdef MOAB_HAVE_CGM 00424 const bool part_surf, const bool ghost, 00425 #else 00426 const bool, const bool, 00427 #endif 00428 const int projection_type, const bool recompute_rcb_box, const bool print_time ) 00429 { 00430 // should only be called in serial 00431 if( mbpc->proc_config().proc_size() != 1 ) 00432 { 00433 std::cout << "ZoltanPartitioner::partition_mesh_and_geometry must be called in serial." << std::endl; 00434 return MB_FAILURE; 00435 } 00436 clock_t t = clock(); 00437 if( NULL != zmethod && strcmp( zmethod, "RR" ) && strcmp( zmethod, "RCB" ) && strcmp( zmethod, "RIB" ) && 00438 strcmp( zmethod, "HSFC" ) && strcmp( zmethod, "Hypergraph" ) && strcmp( zmethod, "PHG" ) && 00439 strcmp( zmethod, "PARMETIS" ) && strcmp( zmethod, "OCTPART" ) ) 00440 { 00441 std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be " 00442 << "RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl; 00443 return MB_FAILURE; 00444 } 00445 00446 bool part_geom = false; 00447 if( 0 == strcmp( zmethod, "RR" ) || 0 == strcmp( zmethod, "RCB" ) || 0 == strcmp( zmethod, "RIB" ) || 00448 0 == strcmp( zmethod, "HSFC" ) ) 00449 part_geom = true; // so no adjacency / edges needed 00450 std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB 00451 std::vector< int > ids; // point ids from MOAB 00452 std::vector< int > adjs, length, parts; 00453 std::vector< double > obj_weights, edge_weights; 00454 Range elems; 00455 #ifdef MOAB_HAVE_CGM 00456 DLIList< RefEntity* > entities; 00457 #endif 00458 00459 // Get a mesh from MOAB and diide it across processors. 00460 00461 ErrorCode result = MB_SUCCESS; 00462 00463 // short-circuit everything if RR partition is requested 00464 if( !strcmp( zmethod, "RR" ) ) 00465 { 00466 if( part_geom_mesh_size < 0. ) 00467 { 00468 // get all elements 00469 result = mbImpl->get_entities_by_dimension( 0, part_dim, elems );RR; 00470 if( elems.empty() ) return MB_FAILURE; 00471 // make a trivial assignment vector 00472 std::vector< int > assign_vec( elems.size() ); 00473 int num_per = elems.size() / nparts; 00474 int extra = elems.size() % nparts; 00475 if( extra ) num_per++; 00476 int nstart = 0; 00477 for( int i = 0; i < nparts; i++ ) 00478 { 00479 if( i == extra ) num_per--; 00480 std::fill( &assign_vec[nstart], &assign_vec[nstart + num_per], i ); 00481 nstart += num_per; 00482 } 00483 00484 result = write_partition( nparts, elems, &assign_vec[0], write_as_sets, write_as_tags );RR; 00485 } 00486 else 00487 { 00488 #ifdef MOAB_HAVE_CGM 00489 result = partition_round_robin( nparts );RR; 00490 #endif 00491 } 00492 00493 return result; 00494 } 00495 00496 std::cout << "Assembling graph..." << std::endl; 00497 00498 if( part_geom_mesh_size < 0. ) 00499 { 00500 // if (!part_geom) { 00501 result = assemble_graph( part_dim, pts, ids, adjs, length, elems, part_geom, projection_type );RR; 00502 } 00503 else 00504 { 00505 #ifdef MOAB_HAVE_CGM 00506 result = assemble_graph( part_dim, pts, ids, adjs, length, obj_weights, edge_weights, parts, entities, 00507 part_geom_mesh_size, nparts );RR; 00508 00509 if( debug ) 00510 { 00511 int n_ids = ids.size(); 00512 entities.reset(); 00513 int i_leng = 0; 00514 for( int i = 0; i < n_ids; i++ ) 00515 { 00516 std::cout << "graph_input_ids[" << i << "]=" << ids[i] << ",obj_weights=" << obj_weights[i] 00517 << ",entity_id=" << entities.get_and_step()->id() << ",part=" << parts[i] << std::endl; 00518 for( int j = 0; j < length[i]; j++ ) 00519 { 00520 std::cout << "adjs[" << j << "]=" << adjs[i_leng] << ",edge_weights=" << edge_weights[i_leng] 00521 << std::endl; 00522 i_leng++; 00523 } 00524 } 00525 } 00526 #endif 00527 } 00528 if( print_time ) 00529 { 00530 std::cout << " time to assemble graph: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n"; 00531 t = clock(); 00532 } 00533 double* o_wgt = NULL; 00534 double* e_wgt = NULL; 00535 if( obj_weights.size() > 0 ) o_wgt = &obj_weights[0]; 00536 if( edge_weights.size() > 0 ) e_wgt = &edge_weights[0]; 00537 00538 myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0], o_wgt, e_wgt, &parts[0], 00539 part_geom ); 00540 00541 // Initialize Zoltan. This is a C call. The simple C++ code 00542 // that creates Zoltan objects does not keep track of whether 00543 // Zoltan_Initialize has been called. 00544 00545 if( print_time ) 00546 { 00547 std::cout << " time to initialize points: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n"; 00548 t = clock(); 00549 } 00550 float version; 00551 00552 std::cout << "Initializing zoltan..." << std::endl; 00553 00554 Zoltan_Initialize( argcArg, argvArg, &version ); 00555 00556 // Create Zoltan object. This calls Zoltan_Create. 00557 if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() ); 00558 00559 if( NULL == zmethod || !strcmp( zmethod, "RCB" ) ) 00560 SetRCB_Parameters(recompute_rcb_box); 00561 else if( !strcmp( zmethod, "RIB" ) ) 00562 SetRIB_Parameters(); 00563 else if( !strcmp( zmethod, "HSFC" ) ) 00564 SetHSFC_Parameters(); 00565 else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) ) 00566 { 00567 if( NULL == other_method || ( other_method[0] == '\0' ) ) 00568 SetHypergraph_Parameters( "auto" ); 00569 else 00570 SetHypergraph_Parameters( other_method ); 00571 00572 if( imbal_tol ) 00573 { 00574 std::ostringstream str; 00575 str << imbal_tol; 00576 myZZ->Set_Param( "IMBALANCE_TOL", str.str().c_str() ); // no debug messages 00577 } 00578 } 00579 else if( !strcmp( zmethod, "PARMETIS" ) ) 00580 { 00581 if( NULL == other_method ) 00582 SetPARMETIS_Parameters( "RepartGDiffusion" ); 00583 else 00584 SetPARMETIS_Parameters( other_method ); 00585 } 00586 else if( !strcmp( zmethod, "OCTPART" ) ) 00587 { 00588 if( NULL == other_method ) 00589 SetOCTPART_Parameters( "2" ); 00590 else 00591 SetOCTPART_Parameters( other_method ); 00592 } 00593 00594 // set # requested partitions 00595 char buff[10]; 00596 sprintf( buff, "%d", nparts ); 00597 int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff ); 00598 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00599 00600 // request only partition assignments 00601 retval = myZZ->Set_Param( "RETURN_LISTS", "PARTITION ASSIGNMENTS" ); 00602 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00603 00604 if( obj_weight > 0 ) 00605 { 00606 std::ostringstream str; 00607 str << obj_weight; 00608 retval = myZZ->Set_Param( "OBJ_WEIGHT_DIM", str.str().c_str() ); 00609 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00610 } 00611 00612 if( edge_weight > 0 ) 00613 { 00614 std::ostringstream str; 00615 str << edge_weight; 00616 retval = myZZ->Set_Param( "EDGE_WEIGHT_DIM", str.str().c_str() ); 00617 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00618 } 00619 00620 // Call backs: 00621 00622 myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL ); 00623 myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL ); 00624 myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL ); 00625 myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL ); 00626 myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL ); 00627 myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL ); 00628 if( part_geom_mesh_size > 0. ) { myZZ->Set_Part_Multi_Fn( mbGetPart, NULL ); } 00629 00630 // Perform the load balancing partitioning 00631 00632 int changes; 00633 int numGidEntries; 00634 int numLidEntries; 00635 int dumnum1; 00636 ZOLTAN_ID_PTR dum_local, dum_global; 00637 int *dum1, *dum2; 00638 int num_assign; 00639 ZOLTAN_ID_PTR assign_gid, assign_lid; 00640 int *assign_procs, *assign_parts; 00641 00642 std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nparts 00643 << " processors..." << std::endl; 00644 #ifndef NDEBUG 00645 #if 0 00646 if (NULL == zmethod || !strcmp(zmethod, "RCB")) 00647 Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 1, 0, 0); 00648 00649 if ( !strcmp(zmethod, "PHG")) 00650 Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 0, 1, 1); 00651 #endif 00652 #endif 00653 retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, dumnum1, dum_global, dum_local, dum1, dum2, 00654 num_assign, assign_gid, assign_lid, assign_procs, assign_parts ); 00655 if( ZOLTAN_OK != retval ) return MB_FAILURE; 00656 00657 if( print_time ) 00658 { 00659 std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n"; 00660 t = clock(); 00661 } 00662 // take results & write onto MOAB partition sets 00663 std::cout << "Saving partition information to MOAB..." << std::endl; 00664 00665 if( part_geom_mesh_size < 0. ) 00666 { 00667 // if (!part_geom) { 00668 result = write_partition( nparts, elems, assign_parts, write_as_sets, write_as_tags ); 00669 } 00670 else 00671 { 00672 #ifdef MOAB_HAVE_CGM 00673 result = write_partition( nparts, entities, assign_parts, obj_weights, part_surf, ghost ); 00674 #endif 00675 } 00676 00677 if( print_time ) 00678 { 00679 std::cout << " time to write partition in memory " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n"; 00680 t = clock(); 00681 } 00682 00683 if( MB_SUCCESS != result ) return result; 00684 00685 // Free the memory allocated for lists returned by LB_Parition() 00686 myZZ->LB_Free_Part( &assign_gid, &assign_lid, &assign_procs, &assign_parts ); 00687 00688 return MB_SUCCESS; 00689 } 00690 00691 ErrorCode ZoltanPartitioner::include_closure() 00692 { 00693 ErrorCode result; 00694 Range ents; 00695 Range adjs; 00696 std::cout << "Adding closure..." << std::endl; 00697 00698 for( Range::iterator rit = partSets.begin(); rit != partSets.end(); ++rit ) 00699 { 00700 00701 // get the top-dimensional entities in the part 00702 result = mbImpl->get_entities_by_handle( *rit, ents, true );RR; 00703 00704 if( ents.empty() ) continue; 00705 00706 // get intermediate-dimensional adjs and add to set 00707 for( int d = mbImpl->dimension_from_handle( *ents.begin() ) - 1; d >= 1; d-- ) 00708 { 00709 adjs.clear(); 00710 result = mbImpl->get_adjacencies( ents, d, false, adjs, Interface::UNION );RR; 00711 result = mbImpl->add_entities( *rit, adjs );RR; 00712 } 00713 00714 // now get vertices and add to set; only need to do for ents, not for adjs 00715 adjs.clear(); 00716 result = mbImpl->get_adjacencies( ents, 0, false, adjs, Interface::UNION );RR; 00717 result = mbImpl->add_entities( *rit, adjs );RR; 00718 00719 ents.clear(); 00720 } 00721 00722 // now go over non-part entity sets, looking for contained entities 00723 Range sets, part_ents; 00724 result = mbImpl->get_entities_by_type( 0, MBENTITYSET, sets );RR; 00725 for( Range::iterator rit = sets.begin(); rit != sets.end(); ++rit ) 00726 { 00727 // skip parts 00728 if( partSets.find( *rit ) != partSets.end() ) continue; 00729 00730 // get entities in this set, recursively 00731 ents.clear(); 00732 result = mbImpl->get_entities_by_handle( *rit, ents, true );RR; 00733 00734 // now check over all parts 00735 for( Range::iterator rit2 = partSets.begin(); rit2 != partSets.end(); ++rit2 ) 00736 { 00737 part_ents.clear(); 00738 result = mbImpl->get_entities_by_handle( *rit2, part_ents, false );RR; 00739 Range int_range = intersect( ents, part_ents ); 00740 if( !int_range.empty() ) 00741 { 00742 // non-empty intersection, add to part set 00743 result = mbImpl->add_entities( *rit2, &( *rit ), 1 );RR; 00744 } 00745 } 00746 } 00747 00748 // finally, mark all the part sets as having the closure 00749 Tag closure_tag; 00750 result = 00751 mbImpl->tag_get_handle( "INCLUDES_CLOSURE", 1, MB_TYPE_INTEGER, closure_tag, MB_TAG_SPARSE | MB_TAG_CREAT );RR; 00752 00753 std::vector< int > closure_vals( partSets.size(), 1 ); 00754 result = mbImpl->tag_set_data( closure_tag, partSets, &closure_vals[0] );RR; 00755 00756 return MB_SUCCESS; 00757 } 00758 00759 ErrorCode ZoltanPartitioner::assemble_graph( const int dimension, std::vector< double >& coords, 00760 std::vector< int >& moab_ids, std::vector< int >& adjacencies, 00761 std::vector< int >& length, Range& elems, bool part_geom, 00762 int projection_type ) 00763 { 00764 // assemble a graph with vertices equal to elements of specified dimension, edges 00765 // signified by list of other elements to which an element is connected 00766 00767 // get the elements of that dimension 00768 ErrorCode result = mbImpl->get_entities_by_dimension( 0, dimension, elems ); 00769 if( MB_SUCCESS != result || elems.empty() ) return result; 00770 00771 // assign global ids 00772 if( assign_global_ids ) 00773 { 00774 EntityHandle rootset = 0; 00775 result = mbpc->assign_global_ids( rootset, dimension, 1, true, true );RR; 00776 } 00777 00778 // now assemble the graph, calling MeshTopoUtil to get bridge adjacencies through d-1 00779 // dimensional neighbors 00780 MeshTopoUtil mtu( mbImpl ); 00781 Range adjs; 00782 // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited 00783 // by MBCN 00784 int neighbors[5 * MAX_SUB_ENTITIES]; 00785 double avg_position[3]; 00786 int moab_id; 00787 00788 // get the global id tag handle 00789 Tag gid = mbImpl->globalId_tag(); 00790 00791 for( Range::iterator rit = elems.begin(); rit != elems.end(); ++rit ) 00792 { 00793 00794 if( !part_geom ) 00795 { 00796 // get bridge adjacencies 00797 adjs.clear(); 00798 result = mtu.get_bridge_adjacencies( *rit, ( dimension > 0 ? dimension - 1 : 3 ), dimension, adjs );RR; 00799 00800 // get the graph vertex ids of those 00801 if( !adjs.empty() ) 00802 { 00803 assert( adjs.size() < 5 * MAX_SUB_ENTITIES ); 00804 result = mbImpl->tag_get_data( gid, adjs, neighbors );RR; 00805 } 00806 00807 // copy those into adjacencies vector 00808 length.push_back( (int)adjs.size() ); 00809 std::copy( neighbors, neighbors + adjs.size(), std::back_inserter( adjacencies ) ); 00810 } 00811 00812 // get average position of vertices 00813 result = mtu.get_average_position( *rit, avg_position );RR; 00814 00815 // get the graph vertex id for this element 00816 result = mbImpl->tag_get_data( gid, &( *rit ), 1, &moab_id );RR; 00817 00818 // copy those into coords vector 00819 moab_ids.push_back( moab_id ); 00820 // transform coordinates to spherical coordinates, if requested 00821 if (projection_type > 0) 00822 IntxUtils::transform_coordinates(avg_position, projection_type); 00823 00824 std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) ); 00825 } 00826 00827 if( debug ) 00828 { 00829 std::cout << "Length vector: " << std::endl; 00830 std::copy( length.begin(), length.end(), std::ostream_iterator< int >( std::cout, ", " ) ); 00831 std::cout << std::endl; 00832 std::cout << "Adjacencies vector: " << std::endl; 00833 std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( std::cout, ", " ) ); 00834 std::cout << std::endl; 00835 std::cout << "Moab_ids vector: " << std::endl; 00836 std::copy( moab_ids.begin(), moab_ids.end(), std::ostream_iterator< int >( std::cout, ", " ) ); 00837 std::cout << std::endl; 00838 std::cout << "Coords vector: " << std::endl; 00839 std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( std::cout, ", " ) ); 00840 std::cout << std::endl; 00841 } 00842 00843 return MB_SUCCESS; 00844 } 00845 00846 #ifdef MOAB_HAVE_CGM 00847 ErrorCode ZoltanPartitioner::assemble_graph( const int /* dimension */, std::vector< double >& /* coords */, 00848 std::vector< int >& moab_ids, std::vector< int >& adjacencies, 00849 std::vector< int >& length, std::vector< double >& obj_weights, 00850 std::vector< double >& edge_weights, std::vector< int >& parts, 00851 DLIList< RefEntity* >& entities, const double part_geom_mesh_size, 00852 const int n_part ) 00853 { 00854 // get body vertex weights 00855 DLIList< RefEntity* > body_list; 00856 gti->ref_entity_list( "body", body_list, CUBIT_FALSE ); 00857 DLIList< RefFace* > all_shared_surfs; 00858 int n_bodies = body_list.size(); 00859 int body_map_index = 1; 00860 int surf_map_index = n_bodies + 1; 00861 int n_bodies_proc = n_bodies / n_part; // # of entities per processor 00862 int i_bodies_proc = n_bodies_proc; // entity index limit for each processor 00863 int proc = 0; 00864 00865 body_list.reset(); 00866 for( int i = 0; i < n_bodies; i++ ) 00867 { 00868 // assign part in round-robin 00869 if( i == i_bodies_proc ) 00870 { 00871 proc++; 00872 if( proc < n_part ) 00873 i_bodies_proc += n_bodies_proc; 00874 else 00875 { 00876 proc %= n_part; 00877 i_bodies_proc++; 00878 } 00879 } 00880 parts.push_back( proc ); 00881 00882 // volume mesh generation load estimation 00883 RefEntity* body = body_list.get_and_step(); 00884 double vertex_w = body->measure(); 00885 vertex_w /= part_geom_mesh_size * part_geom_mesh_size * part_geom_mesh_size; 00886 vertex_w = 3.8559e-5 * vertex_w * log( vertex_w ); 00887 00888 // add child non-interface face weights 00889 DLIList< RefFace* > shared_surfs; 00890 DLIList< RefFace* > faces; 00891 ( dynamic_cast< TopologyEntity* >( body ) )->ref_faces( faces ); 00892 int n_face = faces.size(); 00893 faces.reset(); 00894 for( int j = 0; j < n_face; j++ ) 00895 { 00896 RefFace* face = faces.get_and_step(); 00897 TopologyEntity* te = CAST_TO( face, TopologyEntity ); 00898 if( te->bridge_manager()->number_of_bridges() > 1 ) 00899 { // shared 00900 shared_surfs.append( face ); 00901 } 00902 else 00903 { // non-shared 00904 vertex_w += estimate_face_mesh_load( face, face->measure() ); 00905 } 00906 } 00907 00908 int temp_index = body_map_index++; 00909 body_vertex_map[body->id()] = temp_index; 00910 moab_ids.push_back( temp_index ); // add body as graph vertex 00911 obj_weights.push_back( vertex_w ); // add weight for body graph vertex 00912 00913 if( debug ) 00914 { 00915 std::cout << "body=" << body->id() << ",graph_id=" << temp_index << ",weight=" << vertex_w 00916 << ",part=" << proc << std::endl; 00917 } 00918 00919 // treat shared surfaces connected to this body 00920 int n_shared = shared_surfs.size(); 00921 shared_surfs.reset(); 00922 for( int k = 0; k < n_shared; k++ ) 00923 { // add adjacencies 00924 RefFace* face = shared_surfs.get_and_step(); 00925 std::map< int, int >::iterator iter = surf_vertex_map.find( face->id() ); 00926 if( iter != surf_vertex_map.end() ) { temp_index = ( *iter ).second; } 00927 else 00928 { 00929 temp_index = surf_map_index++; 00930 surf_vertex_map[face->id()] = temp_index; 00931 all_shared_surfs.append( face ); 00932 } 00933 adjacencies.push_back( temp_index ); 00934 double tmp_sw = estimate_face_comm_load( face, part_geom_mesh_size ); 00935 edge_weights.push_back( tmp_sw ); 00936 00937 if( debug ) { std::cout << "adjac=" << temp_index << ",weight=" << tmp_sw << std::endl; } 00938 } 00939 length.push_back( n_shared ); 00940 } 00941 entities += body_list; 00942 00943 // add shared surface as graph vertex 00944 int n_all_shared_surf = all_shared_surfs.size(); 00945 all_shared_surfs.reset(); 00946 for( int i = 0; i < n_all_shared_surf; i++ ) 00947 { 00948 RefEntity* face = all_shared_surfs.get_and_step(); 00949 moab_ids.push_back( surf_vertex_map[face->id()] ); 00950 double face_mesh_load = estimate_face_mesh_load( face, part_geom_mesh_size ); 00951 obj_weights.push_back( face_mesh_load ); 00952 00953 // set surface object graph edges 00954 int min_ind = -1; 00955 double min_load = std::numeric_limits< double >::max(); 00956 ; 00957 DLIList< Body* > parents; 00958 dynamic_cast< TopologyEntity* >( face )->bodies( parents ); 00959 int n_parents = parents.size(); 00960 parents.reset(); 00961 for( int j = 0; j < n_parents; j++ ) 00962 { 00963 int body_gid = body_vertex_map[parents.get_and_step()->id()]; 00964 adjacencies.push_back( body_gid ); // add adjacency 00965 edge_weights.push_back( estimate_face_comm_load( face, part_geom_mesh_size ) ); 00966 00967 if( obj_weights[body_gid - 1] < min_load ) 00968 { 00969 min_ind = body_gid - 1; 00970 min_load = obj_weights[min_ind]; 00971 } 00972 } 00973 length.push_back( n_parents ); 00974 entities.append( dynamic_cast< RefEntity* >( face ) ); 00975 obj_weights[min_ind] += face_mesh_load; 00976 parts.push_back( parts[min_ind] ); 00977 00978 if( debug ) 00979 { 00980 std::cout << "shared_surf=" << face->id() << ",graph_id=" << surf_vertex_map[face->id()] 00981 << ",weight=" << face_mesh_load << ",part=" << parts[min_ind] << std::endl; 00982 } 00983 } 00984 00985 for( size_t i = 0; i < obj_weights.size(); i++ ) 00986 if( obj_weights[i] < 1. ) obj_weights[i] = 1.; 00987 for( size_t i = 0; i < edge_weights.size(); i++ ) 00988 if( edge_weights[i] < 1. ) edge_weights[i] = 1.; 00989 00990 return MB_SUCCESS; 00991 } 00992 00993 double ZoltanPartitioner::estimate_face_mesh_load( RefEntity* face, const double h ) 00994 { 00995 GeometryType type = CAST_TO( face, RefFace )->geometry_type(); 00996 double n = face->measure() / h / h; 00997 double n_logn = n * log( n ); 00998 00999 if( type == PLANE_SURFACE_TYPE ) { return 1.536168737505151e-4 * n_logn; } 01000 else if( type == SPLINE_SURFACE_TYPE ) 01001 { 01002 return 5.910511018383144e-4 * n_logn; 01003 } 01004 else if( type == CONE_SURFACE_TYPE || type == SPHERE_SURFACE_TYPE || type == TORUS_SURFACE_TYPE ) 01005 { 01006 return 2.352511671418708e-4 * n_logn; 01007 } 01008 01009 return 0.0; 01010 } 01011 01012 double ZoltanPartitioner::estimate_face_comm_load( RefEntity* face, const double h ) 01013 { 01014 double peri = 0.; 01015 #if( ( CGM_MAJOR_VERSION == 14 && CGM_MINOR_VERSION > 2 ) || CGM_MAJOR_VERSION >= 15 ) 01016 DLIList< DLIList< RefEdge* > > ref_edge_loops; 01017 #else 01018 DLIList< DLIList< RefEdge* >* > ref_edge_loops; 01019 #endif 01020 CAST_TO( face, RefFace )->ref_edge_loops( ref_edge_loops ); 01021 ref_edge_loops.reset(); 01022 01023 #if( ( CGM_MAJOR_VERSION == 14 && CGM_MINOR_VERSION > 2 ) || CGM_MAJOR_VERSION >= 15 ) 01024 for( int i = 0; i < ref_edge_loops.size(); i++ ) 01025 { 01026 DLIList< RefEdge* > eloop = ref_edge_loops.get_and_step(); 01027 eloop.reset(); 01028 for( int j = 0; j < eloop.size(); j++ ) 01029 { 01030 peri += eloop.get_and_step()->measure(); 01031 } 01032 } 01033 #else 01034 for( int i = 0; i < ref_edge_loops.size(); i++ ) 01035 { 01036 DLIList< RefEdge* >* eloop = ref_edge_loops.get_and_step(); 01037 eloop->reset(); 01038 for( int j = 0; j < eloop->size(); j++ ) 01039 { 01040 peri += eloop->get_and_step()->measure(); 01041 } 01042 } 01043 #endif 01044 01045 // return 104*face->measure()/sqrt(3)/h/h + 56/3*peri/h; 01046 return ( 104 * face->measure() / sqrt( 3 ) / h / h + 56 / 3 * peri / h ) / 700000.; 01047 } 01048 01049 ErrorCode ZoltanPartitioner::write_partition( const int nparts, DLIList< RefEntity* > entities, const int* assignment, 01050 std::vector< double >& obj_weights, const bool part_surf, 01051 const bool ghost ) 01052 { 01053 ErrorCode result; 01054 01055 // actuate CA_BODIES and turn on auto flag for other attributes 01056 CGMApp::instance()->attrib_manager()->register_attrib_type( CA_BODIES, "bodies", "BODIES", CABodies_creator, 01057 CUBIT_TRUE, CUBIT_TRUE, CUBIT_TRUE, CUBIT_TRUE, 01058 CUBIT_TRUE, CUBIT_FALSE ); 01059 CGMApp::instance()->attrib_manager()->auto_flag( CUBIT_TRUE ); 01060 01061 // set partition info to Body at first 01062 int n_entities = entities.size(); 01063 entities.reset(); 01064 for( int i = 0; i < n_entities; i++ ) 01065 { 01066 int proc = assignment[i]; 01067 DLIList< int > shared_procs; 01068 RefEntity* entity = entities.get_and_step(); 01069 01070 if( entity->entity_type_info() == typeid( Body ) ) 01071 { 01072 shared_procs.append( proc ); 01073 TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01074 if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs ); 01075 01076 if( debug ) { std::cout << "body" << entity->id() << "_is_partitioned_to_p" << proc << std::endl; } 01077 01078 // assign to volumes, it should be removed in future 01079 DLIList< RefVolume* > volumes; 01080 ( dynamic_cast< TopologyEntity* >( entity ) )->ref_volumes( volumes ); 01081 int n_vol = volumes.size(); 01082 volumes.reset(); 01083 for( int j = 0; j < n_vol; j++ ) 01084 { 01085 RefEntity* vol = volumes.get_and_step(); 01086 td_par = (TDParallel*)vol->get_TD( &TDParallel::is_parallel ); 01087 if( td_par == NULL ) td_par = new TDParallel( vol, NULL, &shared_procs ); 01088 } 01089 } 01090 } 01091 01092 // set partition info to shared surfaces 01093 entities.reset(); 01094 for( int i = 0; i < n_entities; i++ ) 01095 { 01096 int proc = assignment[i]; 01097 DLIList< int > shared_procs; 01098 RefEntity* entity = entities.get_and_step(); 01099 01100 if( entity->entity_type_info() == typeid( RefFace ) ) 01101 { // surface 01102 DLIList< Body* > parents; 01103 ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parents ); 01104 int n_parents = parents.size(); 01105 if( n_parents != 2 ) 01106 { // check # of parents 01107 std::cerr << "# of parent bodies of interface surface should be 2." << std::endl; 01108 return MB_FAILURE; 01109 } 01110 shared_procs.append( proc ); // local proc 01111 parents.reset(); 01112 for( int j = 0; j < 2; j++ ) 01113 { 01114 RefEntity* parent = parents.get_and_step(); 01115 TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel ); 01116 01117 if( parent_td == NULL ) 01118 { 01119 std::cerr << "parent body should have balanced process." << std::endl; 01120 return MB_FAILURE; 01121 } 01122 int temp_proc = parent_td->get_charge_proc(); 01123 if( temp_proc != proc ) shared_procs.append( temp_proc ); // remote proc 01124 } 01125 01126 if( shared_procs.size() > 1 ) 01127 { // if it is shared by 2 processors 01128 int merge_id = TDUniqueId::get_unique_id( entity ); 01129 TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01130 if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs, NULL, merge_id, 1 ); 01131 01132 if( debug ) 01133 { 01134 std::cout << "surf" << entity->id() << "_is_partitioned_to_p"; 01135 for( int j = 0; j < shared_procs.size(); j++ ) 01136 { 01137 std::cout << "," << shared_procs[j]; 01138 } 01139 std::cout << std::endl; 01140 } 01141 } 01142 } 01143 } 01144 01145 // do non-shared surface partition too 01146 if( part_surf ) 01147 { 01148 result = partition_surface( nparts, entities, assignment, obj_weights );RR; 01149 } 01150 01151 // partition shared edges and vertex 01152 result = partition_child_entities( 1, nparts, part_surf, ghost );RR; 01153 result = partition_child_entities( 0, nparts, part_surf );RR; 01154 01155 if( debug ) 01156 { 01157 entities.reset(); 01158 for( int i = 0; i < n_entities; i++ ) 01159 { 01160 RefEntity* entity = entities.get_and_step(); 01161 if( entity->entity_type_info() == typeid( Body ) ) 01162 { 01163 TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01164 CubitString st = entity->entity_name(); 01165 DLIList< int >* sp = td_par->get_shared_proc_list(); 01166 int n_sp = sp->size(); 01167 sp->reset(); 01168 for( int j = 0; j < n_sp; j++ ) 01169 { 01170 std::cout << "partitioned_" << st.c_str() << ",proc=" << sp->get_and_step() << std::endl; 01171 } 01172 DLIList< int >* gp = td_par->get_ghost_proc_list(); 01173 int n_gp = gp->size(); 01174 sp->reset(); 01175 for( int j = 0; j < n_gp; j++ ) 01176 { 01177 std::cout << "partitioned_" << st.c_str() << ",ghost=" << gp->get_and_step() << std::endl; 01178 } 01179 } 01180 } 01181 } 01182 01183 return MB_SUCCESS; 01184 } 01185 01186 ErrorCode ZoltanPartitioner::partition_surface( const int nparts, DLIList< RefEntity* > entities, const int* assignment, 01187 std::vector< double >& obj_weights ) 01188 { 01189 int i; 01190 double ave_load = 0.; 01191 double* loads = new double[nparts]; 01192 for( i = 0; i < nparts; i++ ) 01193 loads[i] = 0.; 01194 01195 int n_entities = entities.size(); 01196 entities.reset(); 01197 for( i = 0; i < n_entities; i++ ) 01198 { 01199 loads[assignment[i]] += obj_weights[i]; 01200 ave_load += obj_weights[i]; 01201 } 01202 ave_load /= nparts; 01203 01204 if( debug ) 01205 { 01206 for( i = 0; i < nparts; i++ ) 01207 { 01208 std::cout << "loads_before_surface_partition[" << i << "]=" << loads[i] << std::endl; 01209 } 01210 } 01211 01212 int n_iter = 0; 01213 bool b_stop = false; 01214 do 01215 { 01216 // get max and min load processors 01217 int max_proc = nparts, min_proc = 0; 01218 double min_load = std::numeric_limits< double >::max(); 01219 double max_load = std::numeric_limits< double >::min(); 01220 for( i = 0; i < nparts; i++ ) 01221 { 01222 if( loads[i] < min_load ) 01223 { 01224 min_load = loads[i]; 01225 min_proc = i; 01226 } 01227 if( loads[i] > max_load ) 01228 { 01229 max_load = loads[i]; 01230 max_proc = i; 01231 } 01232 } 01233 01234 double load_diff = max_load - ave_load; 01235 if( load_diff > ave_load / 10. ) 01236 { 01237 bool b_moved = false; 01238 entities.reset(); 01239 for( i = 0; i < n_entities; i++ ) 01240 { 01241 if( b_moved ) break; 01242 if( assignment[i] == max_proc && // find maximum load processor bodies 01243 entities[i]->entity_type_info() == typeid( Body ) ) 01244 { 01245 DLIList< RefFace* > faces; 01246 ( dynamic_cast< TopologyEntity* >( entities[i] ) )->ref_faces( faces ); 01247 int n_face = faces.size(); 01248 faces.reset(); 01249 for( int j = 0; j < n_face; j++ ) 01250 { 01251 RefFace* face = faces.get_and_step(); 01252 TopologyEntity* te = CAST_TO( face, TopologyEntity ); 01253 if( te->bridge_manager()->number_of_bridges() < 2 ) 01254 { // non-shared 01255 TDParallel* td_par = (TDParallel*)face->get_TD( &TDParallel::is_parallel ); 01256 if( td_par == NULL ) 01257 { // only consider unpartitioned surfaces 01258 double face_load = face->measure(); 01259 if( load_diff > min_load + face_load - ave_load ) 01260 { 01261 loads[max_proc] -= face_load; 01262 loads[min_proc] += face_load; 01263 int merge_id = TDUniqueId::get_unique_id( face ); 01264 DLIList< int > shared_procs; 01265 shared_procs.append( min_proc ); 01266 shared_procs.append( max_proc ); 01267 td_par = new TDParallel( face, NULL, &shared_procs, NULL, merge_id, 1 ); 01268 01269 if( debug ) 01270 { 01271 std::cout << "non-shared surface " << face->id() << " is moved from p" 01272 << max_proc << " to p" << min_proc << std::endl; 01273 } 01274 b_moved = true; 01275 break; 01276 } 01277 } 01278 } 01279 } 01280 } 01281 } 01282 } 01283 else 01284 b_stop = true; 01285 01286 n_iter++; 01287 } while( !b_stop && n_iter < 50 ); 01288 01289 if( debug ) 01290 { 01291 for( i = 0; i < nparts; i++ ) 01292 { 01293 std::cout << "loads_after_surface_partition[" << i << "]=" << loads[i] << std::endl; 01294 } 01295 } 01296 01297 delete loads; 01298 return MB_SUCCESS; 01299 } 01300 01301 ErrorCode ZoltanPartitioner::partition_round_robin( const int n_part ) 01302 { 01303 int i, j, k; 01304 double* loads = new double[n_part]; // estimated loads for each processor 01305 double* ve_loads = new double[n_part]; // estimated loads for each processor 01306 for( i = 0; i < n_part; i++ ) 01307 { 01308 loads[i] = 0.0; 01309 ve_loads[i] = 0.0; 01310 } 01311 DLIList< RefEntity* > body_entity_list; 01312 gti->ref_entity_list( "body", body_entity_list, CUBIT_FALSE ); 01313 int n_entity = body_entity_list.size(); 01314 int n_entity_proc = n_entity / n_part; // # of entities per processor 01315 int i_entity_proc = n_entity_proc; // entity index limit for each processor 01316 int proc = 0; 01317 RefEntity* entity; 01318 01319 // assign processors to bodies 01320 body_entity_list.reset(); 01321 for( i = 0; i < n_entity; i++ ) 01322 { 01323 if( i == i_entity_proc ) 01324 { 01325 proc++; 01326 if( proc < n_part ) 01327 i_entity_proc += n_entity_proc; 01328 else 01329 { 01330 proc %= n_part; 01331 i_entity_proc++; 01332 } 01333 } 01334 01335 // assign to bodies 01336 entity = body_entity_list.get_and_step(); 01337 DLIList< int > shared_procs; 01338 shared_procs.append( proc ); 01339 TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01340 if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs ); 01341 loads[proc] += entity->measure(); 01342 01343 // assign to volumes, it should be removed in future 01344 DLIList< RefVolume* > volumes; 01345 ( dynamic_cast< TopologyEntity* >( entity ) )->ref_volumes( volumes ); 01346 int n_vol = volumes.size(); 01347 volumes.reset(); 01348 for( j = 0; j < n_vol; j++ ) 01349 { 01350 RefEntity* vol = volumes.get_and_step(); 01351 td_par = (TDParallel*)vol->get_TD( &TDParallel::is_parallel ); 01352 if( td_par == NULL ) td_par = new TDParallel( vol, NULL, &shared_procs ); 01353 } 01354 01355 // add local surface load 01356 DLIList< RefFace* > faces; 01357 ( dynamic_cast< TopologyEntity* >( entity ) )->ref_faces( faces ); 01358 int n_face = faces.size(); 01359 faces.reset(); 01360 for( j = 0; j < n_face; j++ ) 01361 { 01362 RefFace* face = faces.get_and_step(); 01363 TopologyEntity* te = CAST_TO( face, TopologyEntity ); 01364 if( te->bridge_manager()->number_of_bridges() < 2 ) loads[proc] = loads[proc] + face->measure(); 01365 } 01366 01367 // Get all child entities 01368 DLIList< RefEntity* > child_list; 01369 RefEntity::get_all_child_ref_entities( body_entity_list, child_list ); 01370 int n_child = child_list.size(); 01371 01372 // assign processors to interface entities 01373 child_list.reset(); 01374 for( j = 0; j < n_child; j++ ) 01375 { 01376 entity = child_list.get_and_step(); 01377 TopologyEntity* te = CAST_TO( entity, TopologyEntity ); 01378 01379 if( te->bridge_manager()->number_of_bridges() > 1 ) 01380 { 01381 DLIList< Body* > parent_bodies; 01382 DLIList< int > child_shared_procs; // Shared processors of each child entity 01383 ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parent_bodies ); 01384 int n_parent = parent_bodies.size(); 01385 01386 for( k = 0; k < n_parent; k++ ) 01387 { 01388 RefEntity* parent_vol = CAST_TO( parent_bodies.get_and_step(), RefEntity ); 01389 TDParallel* parent_td = (TDParallel*)parent_vol->get_TD( &TDParallel::is_parallel ); 01390 01391 if( parent_td == NULL ) 01392 { 01393 PRINT_ERROR( "parent Volume has to be partitioned." ); 01394 delete[] loads; 01395 delete[] ve_loads; 01396 return MB_FAILURE; 01397 } 01398 child_shared_procs.append_unique( parent_td->get_charge_proc() ); 01399 } 01400 01401 if( child_shared_procs.size() > 1 ) 01402 { // if it is interface 01403 td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01404 if( td_par == NULL ) 01405 { 01406 int merge_id = TDUniqueId::get_unique_id( entity ); 01407 if( entity->entity_type_info() == typeid( RefFace ) ) 01408 { // face 01409 if( child_shared_procs.size() != 2 ) 01410 { 01411 PRINT_ERROR( "Error: # of shared processors of interface surface " 01412 "should be 2." ); 01413 delete[] loads; 01414 delete[] ve_loads; 01415 return MB_FAILURE; 01416 } 01417 01418 // balance interface surface loads 01419 if( loads[child_shared_procs[0]] > loads[child_shared_procs[1]] ) 01420 child_shared_procs.reverse(); 01421 01422 loads[child_shared_procs[0]] = loads[child_shared_procs[0]] + entity->measure(); 01423 td_par = new TDParallel( entity, NULL, &child_shared_procs, NULL, merge_id, 1 ); 01424 } // face 01425 else if( entity->entity_type_info() == typeid( RefEdge ) || 01426 entity->entity_type_info() == typeid( RefVertex ) ) 01427 { // edge or vertex 01428 // balance interface surface loads 01429 int min_p = child_shared_procs[0]; 01430 int n_shared_proc = child_shared_procs.size(); 01431 for( k = 1; k < n_shared_proc; k++ ) 01432 { 01433 if( ve_loads[child_shared_procs[k]] < ve_loads[min_p] ) min_p = child_shared_procs[k]; 01434 } 01435 ve_loads[min_p] = ve_loads[min_p] + entity->measure(); 01436 child_shared_procs.remove( min_p ); 01437 child_shared_procs.insert_first( min_p ); 01438 td_par = new TDParallel( entity, NULL, &child_shared_procs, NULL, merge_id, 1 ); 01439 } // edge or vertex 01440 } // if (td_par == NULL) 01441 } // if it is interface 01442 } // if (te->bridge_manager()->number_of_bridges() > 1) 01443 } // for (j = 0; j < n_child; j++) 01444 } // for (i = 0; i < n_entity; i++) 01445 01446 delete[] loads; 01447 delete[] ve_loads; 01448 01449 return MB_SUCCESS; 01450 } 01451 01452 // partition child entities to one of parent entity shared processors 01453 ErrorCode ZoltanPartitioner::partition_child_entities( const int dim, const int n_part, const bool part_surf, 01454 const bool ghost ) 01455 { 01456 DLIList< RefEntity* > entity_list; 01457 if( dim == 0 ) 01458 gti->ref_entity_list( "vertex", entity_list, CUBIT_FALSE ); 01459 else if( dim == 1 ) 01460 gti->ref_entity_list( "curve", entity_list, CUBIT_FALSE ); 01461 else 01462 { 01463 std::cerr << "Dimention should be from 0 to 1." << std::endl; 01464 return MB_FAILURE; 01465 } 01466 01467 int i, j, k; 01468 int n_entity = entity_list.size(); 01469 double* loads = new double[n_part]; 01470 for( i = 0; i < n_part; i++ ) 01471 loads[i] = 0.; 01472 entity_list.reset(); 01473 01474 for( i = 0; i < n_entity; i++ ) 01475 { // for all entities 01476 RefEntity* entity = entity_list.get_and_step(); 01477 TopologyEntity* te = CAST_TO( entity, TopologyEntity ); 01478 01479 if( !part_surf && te->bridge_manager()->number_of_bridges() < 2 ) continue; 01480 01481 DLIList< int > shared_procs; 01482 DLIList< Body* > parents; 01483 ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parents ); 01484 int n_parents = parents.size(); 01485 std::set< int > s_proc; 01486 parents.reset(); 01487 01488 // get shared processors from parent bodies 01489 for( j = 0; j < n_parents; j++ ) 01490 { 01491 RefEntity* parent = parents.get_and_step(); 01492 TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel ); 01493 01494 if( parent_td != NULL ) 01495 { 01496 DLIList< int >* parent_procs = parent_td->get_shared_proc_list(); 01497 int n_shared = parent_procs->size(); 01498 parent_procs->reset(); 01499 for( k = 0; k < n_shared; k++ ) 01500 { 01501 int p = parent_procs->get_and_step(); 01502 s_proc.insert( p ); 01503 } 01504 } 01505 } 01506 01507 if( part_surf ) 01508 { // also get shared procs from parent surfaces 01509 DLIList< RefFace* > parent_faces; 01510 ( dynamic_cast< TopologyEntity* >( entity ) )->ref_faces( parent_faces ); 01511 int n_pface = parent_faces.size(); 01512 parent_faces.reset(); 01513 01514 // get shared processors from parent faces 01515 for( j = 0; j < n_pface; j++ ) 01516 { 01517 RefEntity* parent = parent_faces.get_and_step(); 01518 TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel ); 01519 01520 if( parent_td != NULL ) 01521 { 01522 DLIList< int >* parent_procs = parent_td->get_shared_proc_list(); 01523 int n_shared = parent_procs->size(); 01524 parent_procs->reset(); 01525 01526 for( k = 0; k < n_shared; k++ ) 01527 { 01528 int p = parent_procs->get_and_step(); 01529 s_proc.insert( p ); 01530 } 01531 } 01532 } 01533 } 01534 01535 // find the minimum load processor and put it 01536 // at the front of the shared_procs list 01537 if( s_proc.size() > 1 ) 01538 { 01539 int min_proc = 0; 01540 double min_load = std::numeric_limits< double >::max(); 01541 std::set< int >::iterator iter = s_proc.begin(); 01542 std::set< int >::iterator end_iter = s_proc.end(); 01543 for( ; iter != end_iter; ++iter ) 01544 { 01545 if( loads[*iter] < min_load ) 01546 { 01547 min_load = loads[*iter]; 01548 min_proc = *iter; 01549 } 01550 } 01551 01552 if( dim == 1 ) 01553 loads[min_proc] += entity->measure(); 01554 else if( dim == 0 ) 01555 loads[min_proc] += 1.; 01556 shared_procs.append( min_proc ); 01557 iter = s_proc.begin(); 01558 end_iter = s_proc.end(); 01559 for( ; iter != end_iter; ++iter ) 01560 { 01561 if( *iter != min_proc ) { shared_procs.append( *iter ); } 01562 } 01563 01564 // add ghost geometries to shared processors for edge 01565 if( ghost ) 01566 { 01567 parents.reset(); 01568 for( j = 0; j < n_parents; j++ ) 01569 { // for all parent bodies 01570 RefEntity* parent_vol = CAST_TO( parents.get_and_step(), RefEntity ); 01571 TDParallel* parent_td = (TDParallel*)parent_vol->get_TD( &TDParallel::is_parallel ); 01572 int n_shared_proc = shared_procs.size(); 01573 01574 for( k = 0; k < n_shared_proc; k++ ) 01575 { 01576 parent_td->add_ghost_proc( shared_procs[k] ); 01577 } 01578 } 01579 } 01580 01581 // initialize tool data 01582 int merge_id = TDUniqueId::get_unique_id( entity ); 01583 TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel ); 01584 if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs, NULL, merge_id, 1 ); 01585 } 01586 } 01587 01588 delete loads; 01589 return MB_SUCCESS; 01590 } 01591 #endif 01592 01593 ErrorCode ZoltanPartitioner::write_partition( const int nparts, Range& elems, const int* assignment, 01594 const bool write_as_sets, const bool write_as_tags ) 01595 { 01596 ErrorCode result; 01597 01598 // get the partition set tag 01599 Tag part_set_tag; 01600 int dum_id = -1, i; 01601 result = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag, 01602 MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );RR; 01603 01604 // get any sets already with this tag, and clear them 01605 Range tagged_sets; 01606 result = 01607 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );RR; 01608 if( !tagged_sets.empty() ) 01609 { 01610 result = mbImpl->clear_meshset( tagged_sets );RR; 01611 if( !write_as_sets ) 01612 { 01613 result = mbImpl->tag_delete_data( part_set_tag, tagged_sets );RR; 01614 } 01615 } 01616 01617 if( write_as_sets ) 01618 { 01619 // first, create partition sets and store in vector 01620 partSets.clear(); 01621 01622 if( nparts > (int)tagged_sets.size() ) 01623 { 01624 // too few partition sets - create missing ones 01625 int num_new = nparts - tagged_sets.size(); 01626 for( i = 0; i < num_new; i++ ) 01627 { 01628 EntityHandle new_set; 01629 result = mbImpl->create_meshset( MESHSET_SET, new_set );RR; 01630 tagged_sets.insert( new_set ); 01631 } 01632 } 01633 else if( nparts < (int)tagged_sets.size() ) 01634 { 01635 // too many partition sets - delete extras 01636 int num_del = tagged_sets.size() - nparts; 01637 for( i = 0; i < num_del; i++ ) 01638 { 01639 EntityHandle old_set = tagged_sets.pop_back(); 01640 result = mbImpl->delete_entities( &old_set, 1 );RR; 01641 } 01642 } 01643 // assign partition sets to vector 01644 partSets.swap( tagged_sets ); 01645 01646 // write a tag to those sets denoting they're partition sets, with a value of the 01647 // proc number 01648 int* dum_ids = new int[nparts]; 01649 for( i = 0; i < nparts; i++ ) 01650 dum_ids[i] = i; 01651 01652 result = mbImpl->tag_set_data( part_set_tag, partSets, dum_ids );RR; 01653 // found out by valgrind when we run mbpart 01654 delete[] dum_ids; 01655 dum_ids = NULL; 01656 01657 // assign entities to the relevant sets 01658 std::vector< EntityHandle > tmp_part_sets; 01659 // int N = (int)elems.size(); 01660 std::copy( partSets.begin(), partSets.end(), std::back_inserter( tmp_part_sets ) ); 01661 /*Range::reverse_iterator riter; 01662 for (i = N-1, riter = elems.rbegin(); riter != elems.rend(); ++riter, i--) { 01663 int assigned_part = assignment[i]; 01664 part_ranges[assigned_part].insert(*riter); 01665 //result = mbImpl->add_entities(tmp_part_sets[assignment[i]], &(*rit), 1); RR; 01666 }*/ 01667 01668 Range::iterator rit; 01669 for( i = 0, rit = elems.begin(); rit != elems.end(); ++rit, i++ ) 01670 { 01671 result = mbImpl->add_entities( tmp_part_sets[assignment[i]], &( *rit ), 1 );RR; 01672 } 01673 /*for (i=0; i<nparts; i++) 01674 { 01675 result = mbImpl->add_entities(tmp_part_sets[i], part_ranges[i]); RR; 01676 } 01677 delete [] part_ranges;*/ 01678 // check for empty sets, warn if there are any 01679 Range empty_sets; 01680 for( rit = partSets.begin(); rit != partSets.end(); ++rit ) 01681 { 01682 int num_ents = 0; 01683 result = mbImpl->get_number_entities_by_handle( *rit, num_ents ); 01684 if( MB_SUCCESS != result || !num_ents ) empty_sets.insert( *rit ); 01685 } 01686 if( !empty_sets.empty() ) 01687 { 01688 std::cout << "WARNING: " << empty_sets.size() << " empty sets in partition: "; 01689 for( rit = empty_sets.begin(); rit != empty_sets.end(); ++rit ) 01690 std::cout << *rit << " "; 01691 std::cout << std::endl; 01692 } 01693 } 01694 01695 if( write_as_tags ) 01696 { 01697 // allocate integer-size partitions 01698 result = mbImpl->tag_set_data( part_set_tag, elems, assignment );RR; 01699 } 01700 01701 return MB_SUCCESS; 01702 } 01703 01704 void ZoltanPartitioner::SetRCB_Parameters(const bool recompute_rcb_box) 01705 { 01706 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nRecursive Coordinate Bisection" << std::endl; 01707 // General parameters: 01708 01709 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01710 myZZ->Set_Param( "LB_METHOD", "RCB" ); // recursive coordinate bisection 01711 // myZZ->Set_Param( "RCB_RECOMPUTE_BOX", "1" ); // recompute RCB box if needed ? 01712 01713 // RCB parameters: 01714 01715 myZZ->Set_Param( "RCB_OUTPUT_LEVEL", "1" ); 01716 myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition so that we can infer partitions 01717 // myZZ->Set_Param("RCB_RECTILINEAR_BLOCKS", "1"); // don't split point on boundary 01718 if (recompute_rcb_box) 01719 myZZ->Set_Param( "RCB_RECOMPUTE_BOX", "1" ); 01720 } 01721 01722 void ZoltanPartitioner::SetRIB_Parameters() 01723 { 01724 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nRecursive Inertial Bisection" << std::endl; 01725 // General parameters: 01726 01727 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01728 myZZ->Set_Param( "LB_METHOD", "RIB" ); // Recursive Inertial Bisection 01729 01730 // RIB parameters: 01731 01732 myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition 01733 myZZ->Set_Param( "AVERAGE_CUTS", "1" ); 01734 } 01735 01736 void ZoltanPartitioner::SetHSFC_Parameters() 01737 { 01738 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nHilbert Space Filling Curve" << std::endl; 01739 // General parameters: 01740 01741 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01742 myZZ->Set_Param( "LB_METHOD", "HSFC" ); // perform Hilbert space filling curve 01743 01744 // HSFC parameters: 01745 01746 myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition 01747 } 01748 01749 void ZoltanPartitioner::SetHypergraph_Parameters( const char* phg_method ) 01750 { 01751 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nHypergraph (or PHG): " << std::endl; 01752 // General parameters: 01753 01754 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01755 myZZ->Set_Param( "LB_METHOD", "Hypergraph" ); // Hypergraph (or PHG) 01756 01757 // Hypergraph (or PHG) parameters: 01758 myZZ->Set_Param( "PHG_COARSEPARTITION_METHOD", phg_method ); // CoarsePartitionMethod 01759 } 01760 01761 void ZoltanPartitioner::SetPARMETIS_Parameters( const char* parmetis_method ) 01762 { 01763 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nPARMETIS: " << parmetis_method << std::endl; 01764 // General parameters: 01765 01766 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01767 myZZ->Set_Param( "LB_METHOD", "PARMETIS" ); // the ParMETIS library 01768 01769 // PARMETIS parameters: 01770 01771 myZZ->Set_Param( "PARMETIS_METHOD", parmetis_method ); // method in the library 01772 } 01773 01774 void ZoltanPartitioner::SetOCTPART_Parameters( const char* oct_method ) 01775 { 01776 if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nOctree Partitioning: " << oct_method << std::endl; 01777 // General parameters: 01778 01779 myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages 01780 myZZ->Set_Param( "LB_METHOD", "OCTPART" ); // octree partitioning 01781 01782 // OCTPART parameters: 01783 01784 myZZ->Set_Param( "OCT_METHOD", oct_method ); // the SFC to be used 01785 myZZ->Set_Param( "OCT_OUTPUT_LEVEL", "3" ); 01786 } 01787 01788 int ZoltanPartitioner::mbInitializePoints( int npts, double* pts, int* ids, int* adjs, int* length, double* obj_weights, 01789 double* edge_weights, int* parts, bool part_geom ) 01790 { 01791 unsigned int i; 01792 int j; 01793 int *numPts, *nborProcs = NULL; 01794 int sum, ptsPerProc, ptsAssigned, mySize; 01795 MPI_Status stat; 01796 double* sendPts; 01797 int* sendIds; 01798 int* sendEdges = NULL; 01799 int* sendNborId = NULL; 01800 int* sendProcs; 01801 01802 if( mbpc->proc_config().proc_rank() == 0 ) 01803 { 01804 /* divide pts to start */ 01805 01806 numPts = (int*)malloc( sizeof( int ) * mbpc->proc_config().proc_size() ); 01807 ptsPerProc = npts / mbpc->proc_config().proc_size(); 01808 ptsAssigned = 0; 01809 01810 for( i = 0; i < mbpc->proc_config().proc_size() - 1; i++ ) 01811 { 01812 numPts[i] = ptsPerProc; 01813 ptsAssigned += ptsPerProc; 01814 } 01815 01816 numPts[mbpc->proc_config().proc_size() - 1] = npts - ptsAssigned; 01817 01818 mySize = numPts[mbpc->proc_config().proc_rank()]; 01819 sendPts = pts + ( 3 * numPts[0] ); 01820 sendIds = ids + numPts[0]; 01821 sum = 0; // possible no adjacency sent 01822 if( !part_geom ) 01823 { 01824 sendEdges = length + numPts[0]; 01825 01826 for( j = 0; j < numPts[0]; j++ ) 01827 sum += length[j]; 01828 01829 sendNborId = adjs + sum; 01830 01831 for( j = numPts[0]; j < npts; j++ ) 01832 sum += length[j]; 01833 01834 nborProcs = (int*)malloc( sizeof( int ) * sum ); 01835 } 01836 for( j = 0; j < sum; j++ ) 01837 if( ( i = adjs[j] / ptsPerProc ) < mbpc->proc_config().proc_size() ) 01838 nborProcs[j] = i; 01839 else 01840 nborProcs[j] = mbpc->proc_config().proc_size() - 1; 01841 01842 sendProcs = nborProcs + ( sendNborId - adjs ); 01843 01844 for( i = 1; i < mbpc->proc_config().proc_size(); i++ ) 01845 { 01846 MPI_Send( &numPts[i], 1, MPI_INT, i, 0x00, mbpc->comm() ); 01847 MPI_Send( sendPts, 3 * numPts[i], MPI_DOUBLE, i, 0x01, mbpc->comm() ); 01848 MPI_Send( sendIds, numPts[i], MPI_INT, i, 0x03, mbpc->comm() ); 01849 MPI_Send( sendEdges, numPts[i], MPI_INT, i, 0x06, mbpc->comm() ); 01850 sum = 0; 01851 01852 for( j = 0; j < numPts[i]; j++ ) 01853 sum += sendEdges[j]; 01854 01855 MPI_Send( sendNborId, sum, MPI_INT, i, 0x07, mbpc->comm() ); 01856 MPI_Send( sendProcs, sum, MPI_INT, i, 0x08, mbpc->comm() ); 01857 sendPts += ( 3 * numPts[i] ); 01858 sendIds += numPts[i]; 01859 sendEdges += numPts[i]; 01860 sendNborId += sum; 01861 sendProcs += sum; 01862 } 01863 01864 free( numPts ); 01865 } 01866 else 01867 { 01868 MPI_Recv( &mySize, 1, MPI_INT, 0, 0x00, mbpc->comm(), &stat ); 01869 pts = (double*)malloc( sizeof( double ) * 3 * mySize ); 01870 ids = (int*)malloc( sizeof( int ) * mySize ); 01871 length = (int*)malloc( sizeof( int ) * mySize ); 01872 if( obj_weights != NULL ) obj_weights = (double*)malloc( sizeof( double ) * mySize ); 01873 MPI_Recv( pts, 3 * mySize, MPI_DOUBLE, 0, 0x01, mbpc->comm(), &stat ); 01874 MPI_Recv( ids, mySize, MPI_INT, 0, 0x03, mbpc->comm(), &stat ); 01875 MPI_Recv( length, mySize, MPI_INT, 0, 0x06, mbpc->comm(), &stat ); 01876 sum = 0; 01877 01878 for( j = 0; j < mySize; j++ ) 01879 sum += length[j]; 01880 01881 adjs = (int*)malloc( sizeof( int ) * sum ); 01882 if( edge_weights != NULL ) edge_weights = (double*)malloc( sizeof( double ) * sum ); 01883 nborProcs = (int*)malloc( sizeof( int ) * sum ); 01884 MPI_Recv( adjs, sum, MPI_INT, 0, 0x07, mbpc->comm(), &stat ); 01885 MPI_Recv( nborProcs, sum, MPI_INT, 0, 0x08, mbpc->comm(), &stat ); 01886 } 01887 01888 Points = pts; 01889 GlobalIds = ids; 01890 NumPoints = mySize; 01891 NumEdges = length; 01892 NborGlobalId = adjs; 01893 NborProcs = nborProcs; 01894 ObjWeights = obj_weights; 01895 EdgeWeights = edge_weights; 01896 Parts = parts; 01897 01898 return mySize; 01899 } 01900 01901 void ZoltanPartitioner::mbFinalizePoints( int npts, int numExport, ZOLTAN_ID_PTR exportLocalIDs, int* exportProcs, 01902 int** assignment ) 01903 { 01904 int* MyAssignment; 01905 int i; 01906 int numPts; 01907 MPI_Status stat; 01908 int* recvA; 01909 01910 /* assign pts to start */ 01911 01912 if( mbpc->proc_config().proc_rank() == 0 ) 01913 MyAssignment = (int*)malloc( sizeof( int ) * npts ); 01914 else 01915 MyAssignment = (int*)malloc( sizeof( int ) * NumPoints ); 01916 01917 for( i = 0; i < NumPoints; i++ ) 01918 MyAssignment[i] = mbpc->proc_config().proc_rank(); 01919 01920 for( i = 0; i < numExport; i++ ) 01921 MyAssignment[exportLocalIDs[i]] = exportProcs[i]; 01922 01923 if( mbpc->proc_config().proc_rank() == 0 ) 01924 { 01925 /* collect pts */ 01926 recvA = MyAssignment + NumPoints; 01927 01928 for( i = 1; i < (int)mbpc->proc_config().proc_size(); i++ ) 01929 { 01930 MPI_Recv( &numPts, 1, MPI_INT, i, 0x04, mbpc->comm(), &stat ); 01931 MPI_Recv( recvA, numPts, MPI_INT, i, 0x05, mbpc->comm(), &stat ); 01932 recvA += numPts; 01933 } 01934 01935 *assignment = MyAssignment; 01936 } 01937 else 01938 { 01939 MPI_Send( &NumPoints, 1, MPI_INT, 0, 0x04, mbpc->comm() ); 01940 MPI_Send( MyAssignment, NumPoints, MPI_INT, 0, 0x05, mbpc->comm() ); 01941 free( MyAssignment ); 01942 } 01943 } 01944 01945 int ZoltanPartitioner::mbGlobalSuccess( int rc ) 01946 { 01947 int fail = 0; 01948 unsigned int i; 01949 int* vals = (int*)malloc( mbpc->proc_config().proc_size() * sizeof( int ) ); 01950 01951 MPI_Allgather( &rc, 1, MPI_INT, vals, 1, MPI_INT, mbpc->comm() ); 01952 01953 for( i = 0; i < mbpc->proc_config().proc_size(); i++ ) 01954 { 01955 if( vals[i] != ZOLTAN_OK ) 01956 { 01957 if( 0 == mbpc->proc_config().proc_rank() ) { mbShowError( vals[i], "Result on process " ); } 01958 fail = 1; 01959 } 01960 } 01961 01962 free( vals ); 01963 return fail; 01964 } 01965 01966 void ZoltanPartitioner::mbPrintGlobalResult( const char* s, int begin, int import, int exp, int change ) 01967 { 01968 unsigned int i; 01969 int* v1 = (int*)malloc( 4 * sizeof( int ) ); 01970 int* v2 = NULL; 01971 int* v; 01972 01973 v1[0] = begin; 01974 v1[1] = import; 01975 v1[2] = exp; 01976 v1[3] = change; 01977 01978 if( mbpc->proc_config().proc_rank() == 0 ) 01979 { v2 = (int*)malloc( 4 * mbpc->proc_config().proc_size() * sizeof( int ) ); } 01980 01981 MPI_Gather( v1, 4, MPI_INT, v2, 4, MPI_INT, 0, mbpc->comm() ); 01982 01983 if( mbpc->proc_config().proc_rank() == 0 ) 01984 { 01985 fprintf( stdout, "======%s======\n", s ); 01986 for( i = 0, v = v2; i < mbpc->proc_config().proc_size(); i++, v += 4 ) 01987 { 01988 fprintf( stdout, "%u: originally had %d, import %d, exp %d, %s\n", i, v[0], v[1], v[2], 01989 v[3] ? "a change of partitioning" : "no change" ); 01990 } 01991 fprintf( stdout, "==========================================\n" ); 01992 fflush( stdout ); 01993 01994 free( v2 ); 01995 } 01996 01997 free( v1 ); 01998 } 01999 02000 void ZoltanPartitioner::mbShowError( int val, const char* s ) 02001 { 02002 if( s ) printf( "%s ", s ); 02003 02004 switch( val ) 02005 { 02006 case ZOLTAN_OK: 02007 printf( "%d: SUCCESSFUL\n", mbpc->proc_config().proc_rank() ); 02008 break; 02009 case ZOLTAN_WARN: 02010 printf( "%d: WARNING\n", mbpc->proc_config().proc_rank() ); 02011 break; 02012 case ZOLTAN_FATAL: 02013 printf( "%d: FATAL ERROR\n", mbpc->proc_config().proc_rank() ); 02014 break; 02015 case ZOLTAN_MEMERR: 02016 printf( "%d: MEMORY ALLOCATION ERROR\n", mbpc->proc_config().proc_rank() ); 02017 break; 02018 default: 02019 printf( "%d: INVALID RETURN CODE\n", mbpc->proc_config().proc_rank() ); 02020 break; 02021 } 02022 } 02023 02024 /********************** 02025 ** call backs 02026 **********************/ 02027 02028 int mbGetNumberOfAssignedObjects( void* /* userDefinedData */, int* err ) 02029 { 02030 *err = 0; 02031 return NumPoints; 02032 } 02033 02034 void mbGetObjectList( void* /* userDefinedData */, int /* numGlobalIds */, int /* numLids */, ZOLTAN_ID_PTR gids, 02035 ZOLTAN_ID_PTR lids, int wgt_dim, float* obj_wgts, int* err ) 02036 { 02037 for( int i = 0; i < NumPoints; i++ ) 02038 { 02039 gids[i] = GlobalIds[i]; 02040 lids[i] = i; 02041 if( wgt_dim > 0 ) obj_wgts[i] = ObjWeights[i]; 02042 } 02043 02044 *err = 0; 02045 } 02046 02047 int mbGetObjectSize( void* /* userDefinedData */, int* err ) 02048 { 02049 *err = 0; 02050 return 3; 02051 } 02052 02053 void mbGetObject( void* /* userDefinedData */, int /* numGlobalIds */, int /* numLids */, int numObjs, 02054 ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR lids, int numDim, double* pts, int* err ) 02055 { 02056 int i, id, id3; 02057 int next = 0; 02058 02059 if( numDim != 3 ) 02060 { 02061 *err = 1; 02062 return; 02063 } 02064 02065 for( i = 0; i < numObjs; i++ ) 02066 { 02067 id = lids[i]; 02068 02069 if( ( id < 0 ) || ( id >= NumPoints ) ) 02070 { 02071 *err = 1; 02072 return; 02073 } 02074 02075 id3 = lids[i] * 3; 02076 02077 pts[next++] = Points[id3]; 02078 pts[next++] = Points[id3 + 1]; 02079 pts[next++] = Points[id3 + 2]; 02080 } 02081 } 02082 02083 void mbGetNumberOfEdges( void* /* userDefinedData */, int /* numGlobalIds */, int /* numLids */, int numObjs, 02084 ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR lids, int* numEdges, int* err ) 02085 { 02086 int i, id; 02087 int next = 0; 02088 02089 for( i = 0; i < numObjs; i++ ) 02090 { 02091 id = lids[i]; 02092 02093 if( ( id < 0 ) || ( id >= NumPoints ) ) 02094 { 02095 *err = 1; 02096 return; 02097 } 02098 02099 numEdges[next++] = NumEdges[id]; 02100 } 02101 } 02102 02103 void mbGetEdgeList( void* /* userDefinedData */, int /* numGlobalIds */, int /* numLids */, int numObjs, 02104 ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR lids, int* /* numEdges */, ZOLTAN_ID_PTR nborGlobalIds, 02105 int* nborProcs, int wgt_dim, float* edge_wgts, int* err ) 02106 { 02107 int i, id, idSum, j; 02108 int next = 0; 02109 02110 for( i = 0; i < numObjs; i++ ) 02111 { 02112 id = lids[i]; 02113 02114 if( ( id < 0 ) || ( id >= NumPoints ) ) 02115 { 02116 *err = 1; 02117 return; 02118 } 02119 02120 idSum = 0; 02121 02122 for( j = 0; j < id; j++ ) 02123 idSum += NumEdges[j]; 02124 02125 for( j = 0; j < NumEdges[id]; j++ ) 02126 { 02127 nborGlobalIds[next] = NborGlobalId[idSum]; 02128 nborProcs[next] = NborProcs[idSum]; 02129 if( wgt_dim > 0 ) edge_wgts[next] = EdgeWeights[idSum]; 02130 next++; 02131 idSum++; 02132 } 02133 } 02134 } 02135 02136 void mbGetPart( void* /* userDefinedData */, int /* numGlobalIds */, int /* numLids */, int numObjs, 02137 ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR lids, int* part, int* err ) 02138 { 02139 int i, id; 02140 int next = 0; 02141 02142 for( i = 0; i < numObjs; i++ ) 02143 { 02144 id = lids[i]; 02145 02146 if( ( id < 0 ) || ( id >= NumPoints ) ) 02147 { 02148 *err = 1; 02149 return; 02150 } 02151 02152 part[next++] = Parts[id]; 02153 } 02154 } 02155 02156 // new methods for partition in parallel, used by migrate in iMOAB 02157 ErrorCode ZoltanPartitioner::partition_owned_cells( Range& primary, ParallelComm* pco, 02158 std::multimap< int, int >& extraGraphEdges, 02159 std::map< int, int > procs, int& numNewPartitions, 02160 std::map< int, Range >& distribution, int met ) 02161 { 02162 // start copy 02163 MeshTopoUtil mtu( mbImpl ); 02164 Range adjs; 02165 02166 std::vector< int > adjacencies; 02167 std::vector< int > ids; 02168 ids.resize( primary.size() ); 02169 std::vector< int > length; 02170 std::vector< double > coords; 02171 std::vector< int > nbor_proc; 02172 // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited 02173 // by MBCN 02174 int neighbors[5 * MAX_SUB_ENTITIES]; 02175 int neib_proc[5 * MAX_SUB_ENTITIES]; 02176 double avg_position[3]; 02177 int moab_id; 02178 int primaryDim = mbImpl->dimension_from_handle( *primary.rbegin() ); 02179 // get the global id tag handle 02180 Tag gid = mbImpl->globalId_tag(); 02181 02182 ErrorCode rval = mbImpl->tag_get_data( gid, primary, &ids[0] );MB_CHK_ERR( rval ); 02183 02184 int rank = pco->rank(); // current rank , will be put on regular neighbors 02185 int i = 0; 02186 for( Range::iterator rit = primary.begin(); rit != primary.end(); ++rit, i++ ) 02187 { 02188 EntityHandle cell = *rit; 02189 // get bridge adjacencies for each cell 02190 if( 1 == met ) 02191 { 02192 adjs.clear(); 02193 rval = mtu.get_bridge_adjacencies( cell, ( primaryDim > 0 ? primaryDim - 1 : 3 ), primaryDim, adjs );MB_CHK_ERR( rval ); 02194 02195 // get the graph vertex ids of those 02196 if( !adjs.empty() ) 02197 { 02198 assert( adjs.size() < 5 * MAX_SUB_ENTITIES ); 02199 rval = mbImpl->tag_get_data( gid, adjs, neighbors );MB_CHK_ERR( rval ); 02200 } 02201 // if adjacent to neighbor partitions, add to the list 02202 int size_adjs = (int)adjs.size(); 02203 moab_id = ids[i]; 02204 for( int k = 0; k < size_adjs; k++ ) 02205 neib_proc[k] = rank; // current rank 02206 if( extraGraphEdges.find( moab_id ) != extraGraphEdges.end() ) 02207 { 02208 // it means that the current cell is adjacent to a cell in another partition ; maybe 02209 // a few 02210 std::pair< std::multimap< int, int >::iterator, std::multimap< int, int >::iterator > ret; 02211 ret = extraGraphEdges.equal_range( moab_id ); 02212 for( std::multimap< int, int >::iterator it = ret.first; it != ret.second; ++it ) 02213 { 02214 int otherID = it->second; 02215 neighbors[size_adjs] = otherID; // the id of the other cell, across partition 02216 neib_proc[size_adjs] = procs[otherID]; // this is how we built this map, the cell id maps to 02217 // what proc it came from 02218 size_adjs++; 02219 } 02220 } 02221 // copy those into adjacencies vector 02222 length.push_back( size_adjs ); 02223 std::copy( neighbors, neighbors + size_adjs, std::back_inserter( adjacencies ) ); 02224 std::copy( neib_proc, neib_proc + size_adjs, std::back_inserter( nbor_proc ) ); 02225 } 02226 else if( 2 == met ) 02227 { 02228 if( TYPE_FROM_HANDLE( cell ) == MBVERTEX ) 02229 { 02230 rval = mbImpl->get_coords( &cell, 1, avg_position );MB_CHK_ERR( rval ); 02231 } 02232 else 02233 { 02234 rval = mtu.get_average_position( cell, avg_position );MB_CHK_ERR( rval ); 02235 } 02236 std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) ); 02237 } 02238 } 02239 02240 #ifdef VERBOSE 02241 std::stringstream ff2; 02242 ff2 << "zoltanInput_" << pco->rank() << ".txt"; 02243 std::ofstream ofs; 02244 ofs.open( ff2.str().c_str(), std::ofstream::out ); 02245 ofs << "Length vector: " << std::endl; 02246 std::copy( length.begin(), length.end(), std::ostream_iterator< int >( ofs, ", " ) ); 02247 ofs << std::endl; 02248 ofs << "Adjacencies vector: " << std::endl; 02249 std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( ofs, ", " ) ); 02250 ofs << std::endl; 02251 ofs << "Moab_ids vector: " << std::endl; 02252 std::copy( ids.begin(), ids.end(), std::ostream_iterator< int >( ofs, ", " ) ); 02253 ofs << std::endl; 02254 ofs << "Coords vector: " << std::endl; 02255 std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( ofs, ", " ) ); 02256 ofs << std::endl; 02257 ofs.close(); 02258 #endif 02259 02260 // these are static var in this file, and used in the callbacks 02261 Points = NULL; 02262 if( 1 != met ) Points = &coords[0]; 02263 GlobalIds = &ids[0]; 02264 NumPoints = (int)ids.size(); 02265 NumEdges = &length[0]; 02266 NborGlobalId = &adjacencies[0]; 02267 NborProcs = &nbor_proc[0]; 02268 ObjWeights = NULL; 02269 EdgeWeights = NULL; 02270 Parts = NULL; 02271 02272 float version; 02273 if( pco->rank() == 0 ) std::cout << "Initializing zoltan..." << std::endl; 02274 02275 Zoltan_Initialize( argcArg, argvArg, &version ); 02276 02277 // Create Zoltan object. This calls Zoltan_Create. 02278 if( NULL == myZZ ) myZZ = new Zoltan( pco->comm() ); 02279 02280 // set # requested partitions 02281 char buff[10]; 02282 sprintf( buff, "%d", numNewPartitions ); 02283 int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff ); 02284 if( ZOLTAN_OK != retval ) return MB_FAILURE; 02285 02286 // request parts assignment 02287 retval = myZZ->Set_Param( "RETURN_LISTS", "PARTS" ); 02288 if( ZOLTAN_OK != retval ) return MB_FAILURE; 02289 02290 myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL ); 02291 myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL ); 02292 // due to a bug in zoltan, if method is graph partitioning, do not pass coordinates!! 02293 if( 2 == met ) 02294 { 02295 myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL ); 02296 myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL ); 02297 SetRCB_Parameters(); // geometry 02298 } 02299 else if( 1 == met ) 02300 { 02301 myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL ); 02302 myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL ); 02303 SetHypergraph_Parameters( "auto" ); 02304 } 02305 02306 // Perform the load balancing partitioning 02307 02308 int changes; 02309 int numGidEntries; 02310 int numLidEntries; 02311 int num_import; 02312 ZOLTAN_ID_PTR import_global_ids, import_local_ids; 02313 int* import_procs; 02314 int* import_to_part; 02315 int num_export; 02316 ZOLTAN_ID_PTR export_global_ids, export_local_ids; 02317 int *assign_procs, *assign_parts; 02318 02319 if( pco->rank() == 0 ) 02320 std::cout << "Computing partition using method (1-graph, 2-geom):" << met << " for " << numNewPartitions 02321 << " parts..." << std::endl; 02322 02323 #ifndef NDEBUG 02324 #if 0 02325 static int counter=0; // it may be possible to call function multiple times in a simulation 02326 // give a way to not overwrite the files 02327 // it should work only with a modified version of Zoltan 02328 std::stringstream basename; 02329 if (1==met) 02330 { 02331 basename << "phg_" << counter++; 02332 Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 0, 1, 0); 02333 } 02334 else if (2==met) 02335 { 02336 basename << "rcb_" << counter++; 02337 Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 1, 0, 0); 02338 } 02339 #endif 02340 #endif 02341 retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids, 02342 import_procs, import_to_part, num_export, export_global_ids, export_local_ids, 02343 assign_procs, assign_parts ); 02344 if( ZOLTAN_OK != retval ) return MB_FAILURE; 02345 02346 #ifdef VERBOSE 02347 std::stringstream ff3; 02348 ff3 << "zoltanOutput_" << pco->rank() << ".txt"; 02349 std::ofstream ofs3; 02350 ofs3.open( ff3.str().c_str(), std::ofstream::out ); 02351 ofs3 << " export elements on rank " << rank << " \n"; 02352 ofs3 << "\t index \t gb_id \t local \t proc \t part \n"; 02353 for( int k = 0; k < num_export; k++ ) 02354 { 02355 ofs3 << "\t" << k << "\t" << export_global_ids[k] << "\t" << export_local_ids[k] << "\t" << assign_procs[k] 02356 << "\t" << assign_parts[k] << "\n"; 02357 } 02358 ofs3.close(); 02359 #endif 02360 02361 // basically each local cell is assigned to a part 02362 02363 assert( num_export == (int)primary.size() ); 02364 for( i = 0; i < num_export; i++ ) 02365 { 02366 EntityHandle cell = primary[export_local_ids[i]]; 02367 distribution[assign_parts[i]].insert( cell ); 02368 } 02369 02370 Zoltan::LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part ); 02371 Zoltan::LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts ); 02372 02373 delete myZZ; 02374 myZZ = NULL; 02375 02376 // clear arrays that were resized locally, to free up local memory 02377 02378 std::vector< int >().swap( adjacencies ); 02379 std::vector< int >().swap( ids ); 02380 std::vector< int >().swap( length ); 02381 std::vector< int >().swap( nbor_proc ); 02382 std::vector< double >().swap( coords ); 02383 02384 return MB_SUCCESS; 02385 }