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