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