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