moab
MBZoltan.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines