MOAB: Mesh Oriented datABase  (version 5.4.1)
MetisPartitioner.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 // Contributed by Lorenzo Alessio Botti (SpaFEDTe)
00017 // This implementation is mostly borrowed from the mbzoltan MOAB partitioning tool
00018 
00019 #include <iostream>
00020 #include <cassert>
00021 #include <sstream>
00022 #include <map>
00023 #include <ctime>
00024 
00025 #include "moab/MetisPartitioner.hpp"
00026 #include "moab/Interface.hpp"
00027 #include "Internals.hpp"
00028 #include "moab/Range.hpp"
00029 #include "moab/WriteUtilIface.hpp"
00030 #include "moab/MeshTopoUtil.hpp"
00031 #include "moab/Skinner.hpp"
00032 #include "MBTagConventions.hpp"
00033 #include "moab/CN.hpp"
00034 
00035 using namespace moab;
00036 
00037 const bool debug = false;
00038 
00039 MetisPartitioner::MetisPartitioner( Interface* impl, const bool use_coords )
00040     : PartitionerBase< idx_t >( impl, use_coords )
00041 
00042 {
00043 }
00044 
00045 MetisPartitioner::~MetisPartitioner() {}
00046 
00047 ErrorCode MetisPartitioner::partition_mesh( const idx_t nparts,
00048                                             const char* method,
00049                                             const int part_dim,
00050                                             const bool write_as_sets,
00051                                             const bool write_as_tags,
00052                                             const bool partition_tagged_sets,
00053                                             const bool partition_tagged_ents,
00054                                             const char* aggregating_tag,
00055                                             const bool print_time )
00056 {
00057 #ifdef MOAB_HAVE_MPI
00058     // should only be called in serial
00059     if( mbpc->proc_config().proc_size() != 1 )
00060     {
00061         std::cout << "MetisPartitioner::partition_mesh_and_geometry must be called in serial." << std::endl;
00062         return MB_FAILURE;
00063     }
00064 #endif
00065 
00066     if( NULL != method && strcmp( method, "ML_RB" ) != 0 && strcmp( method, "ML_KWAY" ) != 0 )
00067     {
00068         std::cout << "ERROR: Method must be "
00069                   << "ML_RB or ML_KWAY" << std::endl;
00070         return MB_FAILURE;
00071     }
00072 
00073     std::vector< double > pts;  // x[0], y[0], z[0], ... from MOAB
00074     std::vector< idx_t > ids;   // poidx_t ids from MOAB
00075     std::vector< idx_t > adjs, parts;
00076     std::vector< idx_t > length;
00077     Range elems;
00078     // Get a mesh from MOAB and diide it across processors.
00079 
00080     clock_t t = clock();
00081 
00082     ErrorCode result;
00083     if( !partition_tagged_sets && !partition_tagged_ents )
00084     {
00085         result = assemble_graph( part_dim, pts, ids, adjs, length, elems );MB_CHK_ERR( result );
00086     }
00087     else if( partition_tagged_sets )
00088     {
00089         result = assemble_taggedsets_graph( part_dim, pts, ids, adjs, length, elems, &( *aggregating_tag ) );MB_CHK_ERR( result );
00090     }
00091     else if( partition_tagged_ents )
00092     {
00093         result = assemble_taggedents_graph( part_dim, pts, ids, adjs, length, elems, &( *aggregating_tag ) );MB_CHK_ERR( result );
00094     }
00095     else
00096     {
00097         MB_SET_ERR( MB_FAILURE, "Either partition tags or sets for Metis partitoner" );
00098     }
00099 
00100     if( print_time )
00101     {
00102         std::cout << " time to assemble graph: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
00103         t = clock();
00104     }
00105 
00106     std::cout << "Computing partition using " << method << " method for " << nparts << " processors..." << std::endl;
00107 
00108     idx_t nelems = length.size() - 1;
00109     idx_t* assign_parts;
00110     assign_parts         = (idx_t*)malloc( sizeof( idx_t ) * nelems );
00111     idx_t nconstraidx_ts = 1;
00112     idx_t edgeCut        = 0;
00113     idx_t nOfPartitions  = static_cast< idx_t >( nparts );
00114     idx_t metis_RESULT;
00115 
00116     if( strcmp( method, "ML_KWAY" ) == 0 )
00117     {
00118         idx_t options[METIS_NOPTIONS];
00119         METIS_SetDefaultOptions( options );
00120         options[METIS_OPTION_CONTIG] = 1;
00121         metis_RESULT = METIS_PartGraphKway( &nelems, &nconstraidx_ts, &length[0], &adjs[0], NULL, NULL, NULL,
00122                                             &nOfPartitions, NULL, NULL, options, &edgeCut, assign_parts );
00123     }
00124     else if( strcmp( method, "ML_RB" ) == 0 )
00125     {
00126         idx_t options[METIS_NOPTIONS];
00127         METIS_SetDefaultOptions( options );
00128         options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;  // CUT
00129         options[METIS_OPTION_IPTYPE]  = METIS_IPTYPE_GROW;  // GROW or RANDOM
00130         options[METIS_OPTION_CTYPE]   = METIS_CTYPE_RM;     // RM or SHEM
00131         options[METIS_OPTION_RTYPE]   = METIS_RTYPE_FM;     // FM
00132         options[METIS_OPTION_NCUTS]   = 10;                 // Number of different partitionings to compute, then
00133                                                             // chooses the best one, default = 1
00134         options[METIS_OPTION_NITER]   = 10;                 // Number of refinements steps, default = 10
00135         options[METIS_OPTION_UFACTOR] = 30;                 // Imabalance, default = 1
00136         options[METIS_OPTION_DBGLVL]  = METIS_DBG_INFO;
00137         metis_RESULT = METIS_PartGraphRecursive( &nelems, &nconstraidx_ts, &length[0], &adjs[0], NULL, NULL, NULL,
00138                                                  &nOfPartitions, NULL, NULL, options, &edgeCut, assign_parts );
00139     }
00140     else
00141         MB_SET_ERR( MB_FAILURE, "Either ML_KWAY or ML_RB needs to be specified for Metis partitioner" );
00142 
00143     if( print_time )
00144     {
00145         std::cout << " time to partition: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
00146         t = clock();
00147     }
00148 
00149 #ifdef MOAB_HAVE_MPI
00150     // assign global node ids, starting from one! TODO
00151     if( assign_global_ids )
00152     {
00153         EntityHandle rootset = 0;
00154         result               = mbpc->assign_global_ids( rootset, part_dim, 1, true, false );MB_CHK_ERR( result );
00155     }
00156 #endif
00157 
00158     if( metis_RESULT != METIS_OK ) return MB_FAILURE;
00159 
00160     // take results & write onto MOAB partition sets
00161     std::cout << "Saving partition information to MOAB..." << std::endl;
00162     {
00163         if( partition_tagged_sets || partition_tagged_ents )
00164         {
00165             result = write_aggregationtag_partition( nparts, elems, assign_parts, write_as_sets, write_as_tags );MB_CHK_ERR( result );
00166         }
00167         else
00168         {
00169             result = write_partition( nparts, elems, assign_parts, write_as_sets, write_as_tags );MB_CHK_ERR( result );
00170         }
00171     }
00172 
00173     if( print_time )
00174     {
00175         std::cout << " time to write partition in memory " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
00176         t = clock();
00177     }
00178     free( assign_parts );
00179 
00180     return MB_SUCCESS;
00181 }
00182 
00183 ErrorCode MetisPartitioner::assemble_taggedents_graph( const int dimension,
00184                                                        std::vector< double >& coords,
00185                                                        std::vector< idx_t >& moab_ids,
00186                                                        std::vector< idx_t >& adjacencies,
00187                                                        std::vector< idx_t >& length,
00188                                                        Range& elems,
00189                                                        const char* aggregating_tag )
00190 {
00191     Tag partSetTag;
00192     ErrorCode result = mbImpl->tag_get_handle( aggregating_tag, 1, MB_TYPE_INTEGER, partSetTag );
00193     if( MB_SUCCESS != result ) return result;
00194 
00195     Range allSubElems;
00196     result = mbImpl->get_entities_by_dimension( 0, dimension, allSubElems );
00197     if( MB_SUCCESS != result || allSubElems.empty() ) return result;
00198     idx_t partSet;
00199     std::map< idx_t, Range > aggloElems;
00200     for( Range::iterator rit = allSubElems.begin(); rit != allSubElems.end(); rit++ )
00201     {
00202         EntityHandle entity = *rit;
00203         result              = mbImpl->tag_get_data( partSetTag, &entity, 1, &partSet );
00204         if( MB_SUCCESS != result ) return result;
00205         if( partSet >= 0 ) aggloElems[partSet].insert( entity );
00206     }
00207     // clear aggregating tag data
00208     TagType type;
00209     result = mbImpl->tag_get_type( partSetTag, type );
00210     if( type == MB_TAG_DENSE )
00211     {
00212         // clear tag on ents and sets
00213         result = mbImpl->tag_delete( partSetTag );
00214         if( MB_SUCCESS != result ) return result;
00215     }
00216     if( type == MB_TAG_SPARSE )
00217     {
00218         // clear tag on ents
00219         result = mbImpl->tag_delete_data( partSetTag, allSubElems );
00220         if( MB_SUCCESS != result ) return result;
00221         // clear tag on sets
00222         result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &partSetTag, 0, 1, elems );
00223         if( MB_SUCCESS != result ) return result;
00224         result = mbImpl->tag_delete_data( partSetTag, elems );
00225         if( MB_SUCCESS != result ) return result;
00226         elems.clear();
00227     }
00228     result =
00229         mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, partSetTag, MB_TAG_SPARSE | MB_TAG_CREAT );
00230     if( MB_SUCCESS != result ) return result;
00231 
00232     for( std::map< idx_t, Range >::iterator mit = aggloElems.begin(); mit != aggloElems.end(); mit++ )
00233     {
00234         EntityHandle new_set;
00235         result = mbImpl->create_meshset( MESHSET_SET, new_set );
00236         if( MB_SUCCESS != result ) return result;
00237         result = mbImpl->add_entities( new_set, mit->second );
00238         if( MB_SUCCESS != result ) return result;
00239         result = mbImpl->tag_set_data( partSetTag, &new_set, 1, &mit->first );
00240         if( MB_SUCCESS != result ) return result;
00241     }
00242 
00243     result =
00244         assemble_taggedsets_graph( dimension, coords, moab_ids, adjacencies, length, elems, &( *aggregating_tag ) );
00245     return MB_SUCCESS;
00246 }
00247 
00248 ErrorCode MetisPartitioner::assemble_taggedsets_graph( const int dimension,
00249                                                        std::vector< double >& coords,
00250                                                        std::vector< idx_t >& moab_ids,
00251                                                        std::vector< idx_t >& adjacencies,
00252                                                        std::vector< idx_t >& length,
00253                                                        Range& elems,
00254                                                        const char* aggregating_tag )
00255 {
00256     length.push_back( 0 );
00257     // assemble a graph with vertices equal to elements of specified dimension, edges
00258     // signified by list of other elements to which an element is connected
00259 
00260     // get the tagged elements
00261     Tag partSetTag;
00262     ErrorCode result = mbImpl->tag_get_handle( aggregating_tag, 1, MB_TYPE_INTEGER, partSetTag );MB_CHK_ERR( result );
00263     // ErrorCode result = mbImpl->tag_get_handle("PARALLEL_PARTITION_SET", 1, MB_TYPE_INTEGER,
00264     // partSetTag);MB_CHK_ERR(result);
00265 
00266     result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &partSetTag, 0, 1, elems );
00267     if( MB_SUCCESS != result || elems.empty() ) return result;
00268 
00269     // assign globla ids to elem sets based on aggregating_tag data
00270     Tag gid_tag;
00271     idx_t zero1 = -1;
00272     result =
00273         mbImpl->tag_get_handle( "GLOBAL_ID_AGGLO", 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_SPARSE | MB_TAG_CREAT, &zero1 );MB_CHK_ERR( result );
00274     for( Range::iterator rit = elems.begin(); rit != elems.end(); rit++ )
00275     {
00276         idx_t partSet;
00277         result = mbImpl->tag_get_data( partSetTag, &( *rit ), 1, &partSet );MB_CHK_ERR( result );
00278         result = mbImpl->tag_set_data( gid_tag, &( *rit ), 1, &partSet );MB_CHK_ERR( result );
00279     }
00280     // clear aggregating tag data
00281     TagType type;
00282     result = mbImpl->tag_get_type( partSetTag, type );MB_CHK_ERR( result );
00283     if( type == MB_TAG_DENSE )
00284     {
00285         result = mbImpl->tag_delete( partSetTag );MB_CHK_ERR( result );
00286     }
00287     if( type == MB_TAG_SPARSE )
00288     {
00289         result = mbImpl->tag_delete_data( partSetTag, elems );MB_CHK_ERR( result );
00290     }
00291 
00292     // assemble the graph, using Skinner to get d-1 dimensional neighbors and then idx_tersecting to
00293     // get adjacencies
00294     std::vector< Range > skin_subFaces( elems.size() );
00295     unsigned int i = 0;
00296     for( Range::iterator rit = elems.begin(); rit != elems.end(); rit++ )
00297     {
00298         Range part_ents;
00299         result = mbImpl->get_entities_by_handle( *rit, part_ents, false );
00300         if( mbImpl->dimension_from_handle( *part_ents.rbegin() ) !=
00301             mbImpl->dimension_from_handle( *part_ents.begin() ) )
00302         {
00303             Range::iterator lower = part_ents.lower_bound( CN::TypeDimensionMap[0].first ),
00304                             upper = part_ents.upper_bound( CN::TypeDimensionMap[dimension - 1].second );
00305             part_ents.erase( lower, upper );
00306         }
00307         Skinner skinner( mbImpl );
00308         result = skinner.find_skin( 0, part_ents, false, skin_subFaces[i], NULL, false, true, false );MB_CHK_ERR( result );
00309         i++;
00310     }
00311     std::vector< EntityHandle > adjs;
00312     std::vector< idx_t > neighbors;
00313     double avg_position[3];
00314     idx_t moab_id;
00315     MeshTopoUtil mtu( mbImpl );
00316     for( unsigned int k = 0; k < i; k++ )
00317     {
00318         // get bridge adjacencies for element k
00319         adjs.clear();
00320         for( unsigned int t = 0; t < i; t++ )
00321         {
00322             if( t != k )
00323             {
00324                 Range subFaces = intersect( skin_subFaces[k], skin_subFaces[t] );
00325                 if( subFaces.size() > 0 ) adjs.push_back( elems[t] );
00326             }
00327         }
00328         if( !adjs.empty() )
00329         {
00330             neighbors.resize( adjs.size() );
00331             result = mbImpl->tag_get_data( gid_tag, &adjs[0], adjs.size(), &neighbors[0] );MB_CHK_ERR( result );
00332         }
00333         // copy those idx_to adjacencies vector
00334         length.push_back( length.back() + (idx_t)adjs.size() );
00335         std::copy( neighbors.begin(), neighbors.end(), std::back_inserter( adjacencies ) );
00336         // get the graph vertex id for this element
00337         const EntityHandle& setk = elems[k];
00338         result                   = mbImpl->tag_get_data( gid_tag, &setk, 1, &moab_id );
00339         moab_ids.push_back( moab_id );
00340         // get average position of vertices
00341         Range part_ents;
00342         result = mbImpl->get_entities_by_handle( elems[k], part_ents, false );MB_CHK_ERR( result );
00343         result = mtu.get_average_position( part_ents, avg_position );MB_CHK_ERR( result );
00344         std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
00345     }
00346     for( unsigned int k = 0; k < i; k++ )
00347     {
00348         for( unsigned int t = 0; t < k; t++ )
00349         {
00350             Range subFaces = intersect( skin_subFaces[k], skin_subFaces[t] );
00351             if( subFaces.size() > 0 ) mbImpl->delete_entities( subFaces );
00352         }
00353     }
00354 
00355     if( debug )
00356     {
00357         std::cout << "Length vector: " << std::endl;
00358         std::copy( length.begin(), length.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00359         std::cout << std::endl;
00360         std::cout << "Adjacencies vector: " << std::endl;
00361         std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00362         std::cout << std::endl;
00363         std::cout << "Moab_ids vector: " << std::endl;
00364         std::copy( moab_ids.begin(), moab_ids.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00365         std::cout << std::endl;
00366         std::cout << "Coords vector: " << std::endl;
00367         std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( std::cout, ", " ) );
00368         std::cout << std::endl;
00369     }
00370     return MB_SUCCESS;
00371 }
00372 
00373 ErrorCode MetisPartitioner::assemble_graph( const int dimension,
00374                                             std::vector< double >& coords,
00375                                             std::vector< idx_t >& moab_ids,
00376                                             std::vector< idx_t >& adjacencies,
00377                                             std::vector< idx_t >& length,
00378                                             Range& elems )
00379 {
00380     length.push_back( 0 );
00381     // assemble a graph with vertices equal to elements of specified dimension, edges
00382     // signified by list of other elements to which an element is connected
00383 
00384     // get the elements of that dimension
00385     ErrorCode result = mbImpl->get_entities_by_dimension( 0, dimension, elems );
00386     if( MB_SUCCESS != result || elems.empty() ) return result;
00387 
00388 #ifdef MOAB_HAVE_MPI
00389     // assign global ids
00390     if( assign_global_ids )
00391     {
00392         result = mbpc->assign_global_ids( 0, dimension, 0 );MB_CHK_ERR( result );
00393     }
00394 #endif
00395 
00396     // now assemble the graph, calling MeshTopoUtil to get bridge adjacencies through d-1
00397     // dimensional neighbors
00398     MeshTopoUtil mtu( mbImpl );
00399     Range adjs;
00400     // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited
00401     // by MBCN
00402     int neighbors[5 * MAX_SUB_ENTITIES];  // these will be now indices in the elems range
00403 
00404     double avg_position[3];
00405     int index_in_elems = 0;
00406 
00407     for( Range::iterator rit = elems.begin(); rit != elems.end(); rit++, index_in_elems++ )
00408     {
00409 
00410         // get bridge adjacencies
00411         adjs.clear();
00412         result = mtu.get_bridge_adjacencies( *rit, ( dimension > 0 ? dimension - 1 : 3 ), dimension, adjs );MB_CHK_ERR( result );
00413 
00414         // get the indices in elems range of those
00415         if( !adjs.empty() )
00416         {
00417             int i = 0;
00418             assert( adjs.size() < 5 * MAX_SUB_ENTITIES );
00419             for( Range::iterator ait = adjs.begin(); ait != adjs.end(); ait++, i++ )
00420             {
00421                 EntityHandle adjEnt = *ait;
00422                 neighbors[i]        = elems.index( adjEnt );
00423             }
00424         }
00425 
00426         // copy those idx_to adjacencies vector
00427         length.push_back( length.back() + (idx_t)adjs.size() );
00428         // conversion made to idx_t
00429         std::copy( neighbors, neighbors + adjs.size(), std::back_inserter( adjacencies ) );
00430 
00431         // get average position of vertices
00432         result = mtu.get_average_position( *rit, avg_position );MB_CHK_ERR( result );
00433 
00434         // get the graph vertex id for this element; it is now index in elems
00435         moab_ids.push_back( index_in_elems );  // conversion made to idx_t
00436 
00437         // copy_to coords vector
00438         std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
00439     }
00440 
00441     if( debug )
00442     {
00443         std::cout << "Length vector: " << std::endl;
00444         std::copy( length.begin(), length.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00445         std::cout << std::endl;
00446         std::cout << "Adjacencies vector: " << std::endl;
00447         std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00448         std::cout << std::endl;
00449         std::cout << "Moab_ids vector: " << std::endl;
00450         std::copy( moab_ids.begin(), moab_ids.end(), std::ostream_iterator< idx_t >( std::cout, ", " ) );
00451         std::cout << std::endl;
00452         std::cout << "Coords vector: " << std::endl;
00453         std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( std::cout, ", " ) );
00454         std::cout << std::endl;
00455     }
00456 
00457     return MB_SUCCESS;
00458 }
00459 
00460 ErrorCode MetisPartitioner::write_aggregationtag_partition( const idx_t nparts,
00461                                                             Range& elems,
00462                                                             const idx_t* assignment,
00463                                                             const bool write_as_sets,
00464                                                             const bool write_as_tags )
00465 {
00466     ErrorCode result;
00467 
00468     // get the partition set tag
00469     Tag part_set_tag;
00470     result =
00471         mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_ERR( result );
00472 
00473     // get any sets already with this tag, and clear them
00474     Range tagged_sets;
00475     result =
00476         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );MB_CHK_ERR( result );
00477     if( !tagged_sets.empty() )
00478     {
00479         result = mbImpl->clear_meshset( tagged_sets );
00480         if( !write_as_sets )
00481         {
00482             result = mbImpl->tag_delete_data( part_set_tag, tagged_sets );MB_CHK_ERR( result );
00483         }
00484     }
00485 
00486     if( write_as_sets )
00487     {
00488         // first, create partition sets and store in vector
00489         partSets.clear();
00490 
00491         if( nparts > (idx_t)tagged_sets.size() )
00492         {
00493             // too few partition sets - create missing ones
00494             idx_t num_new = nparts - tagged_sets.size();
00495             for( idx_t i = 0; i < num_new; i++ )
00496             {
00497                 EntityHandle new_set;
00498                 result = mbImpl->create_meshset( MESHSET_SET, new_set );MB_CHK_ERR( result );
00499                 tagged_sets.insert( new_set );
00500             }
00501         }
00502         else if( nparts < (idx_t)tagged_sets.size() )
00503         {
00504             // too many partition sets - delete extras
00505             idx_t num_del = tagged_sets.size() - nparts;
00506             for( idx_t i = 0; i < num_del; i++ )
00507             {
00508                 EntityHandle old_set = tagged_sets.pop_back();
00509                 result               = mbImpl->delete_entities( &old_set, 1 );MB_CHK_ERR( result );
00510             }
00511         }
00512 
00513         // assign partition sets to vector
00514         partSets.swap( tagged_sets );
00515 
00516         // write a tag to those sets denoting they're partition sets, with a value of the
00517         // proc number
00518         idx_t* dum_ids = new idx_t[nparts];
00519         for( idx_t i = 0; i < nparts; i++ )
00520             dum_ids[i] = i;
00521 
00522         result = mbImpl->tag_set_data( part_set_tag, partSets, dum_ids );MB_CHK_ERR( result );
00523 
00524         // assign entities to the relevant sets
00525         std::vector< EntityHandle > tmp_part_sets;
00526         std::copy( partSets.begin(), partSets.end(), std::back_inserter( tmp_part_sets ) );
00527         Range::iterator rit;
00528         unsigned j = 0;
00529         for( rit = elems.begin(); rit != elems.end(); rit++, j++ )
00530         {
00531             result = mbImpl->add_entities( tmp_part_sets[assignment[j]], &( *rit ), 1 );MB_CHK_ERR( result );
00532         }
00533 
00534         // check for empty sets, warn if there are any
00535         Range empty_sets;
00536         for( rit = partSets.begin(); rit != partSets.end(); rit++ )
00537         {
00538             int num_ents = 0;
00539             result       = mbImpl->get_number_entities_by_handle( *rit, num_ents );
00540             if( MB_SUCCESS != result || !num_ents ) empty_sets.insert( *rit );
00541         }
00542         if( !empty_sets.empty() )
00543         {
00544             std::cout << "WARNING: " << empty_sets.size() << " empty sets in partition: ";
00545             for( rit = empty_sets.begin(); rit != empty_sets.end(); rit++ )
00546                 std::cout << *rit << " ";
00547             std::cout << std::endl;
00548         }
00549     }
00550 
00551     if( write_as_tags )
00552     {
00553         Tag gid_tag;
00554         result = mbImpl->tag_get_handle( "GLOBAL_ID_AGGLO", 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_SPARSE );MB_CHK_ERR( result );
00555 
00556         // allocate idx_teger-size partitions
00557         unsigned int i = 0;
00558         idx_t gid;
00559         for( Range::iterator rit = elems.begin(); rit != elems.end(); rit++ )
00560         {
00561             result = mbImpl->tag_get_data( gid_tag, &( *rit ), 1, &gid );
00562             Range part_ents;
00563             //      std::cout<<"part ents "<<part_ents.size()<<std::endl;
00564             result = mbImpl->get_entities_by_handle( *rit, part_ents, false );MB_CHK_ERR( result );
00565 
00566             for( Range::iterator eit = part_ents.begin(); eit != part_ents.end(); eit++ )
00567             {
00568                 result = mbImpl->tag_set_data( part_set_tag, &( *eit ), 1, &assignment[i] );MB_CHK_ERR( result );
00569 
00570                 result = mbImpl->tag_set_data( gid_tag, &( *eit ), 1, &gid );MB_CHK_ERR( result );
00571             }
00572             i++;
00573         }
00574     }
00575     return MB_SUCCESS;
00576 }
00577 
00578 ErrorCode MetisPartitioner::write_partition( const idx_t nparts,
00579                                              Range& elems,
00580                                              const idx_t* assignment,
00581                                              const bool write_as_sets,
00582                                              const bool write_as_tags )
00583 {
00584     ErrorCode result;
00585 
00586     // get the partition set tag
00587     Tag part_set_tag;
00588     idx_t dum_id = -1, i;
00589     result       = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag,
00590                                            MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );MB_CHK_ERR( result );
00591 
00592     // get any sets already with this tag, and clear them
00593     Range tagged_sets;
00594     result =
00595         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );MB_CHK_ERR( result );
00596     if( !tagged_sets.empty() )
00597     {
00598         result = mbImpl->clear_meshset( tagged_sets );
00599         if( !write_as_sets )
00600         {
00601             result = mbImpl->tag_delete_data( part_set_tag, tagged_sets );MB_CHK_ERR( result );
00602         }
00603     }
00604 
00605     if( write_as_sets )
00606     {
00607         // first, create partition sets and store in vector
00608         partSets.clear();
00609 
00610         if( nparts > (int)tagged_sets.size() )
00611         {
00612             // too few partition sets - create missing ones
00613             idx_t num_new = nparts - tagged_sets.size();
00614             for( i = 0; i < num_new; i++ )
00615             {
00616                 EntityHandle new_set;
00617                 result = mbImpl->create_meshset( MESHSET_SET, new_set );MB_CHK_ERR( result );
00618                 tagged_sets.insert( new_set );
00619             }
00620         }
00621         else if( nparts < (idx_t)tagged_sets.size() )
00622         {
00623             // too many partition sets - delete extras
00624             idx_t num_del = tagged_sets.size() - nparts;
00625             for( i = 0; i < num_del; i++ )
00626             {
00627                 EntityHandle old_set = tagged_sets.pop_back();
00628                 result               = mbImpl->delete_entities( &old_set, 1 );MB_CHK_ERR( result );
00629             }
00630         }
00631 
00632         // assign partition sets to vector
00633         partSets.swap( tagged_sets );
00634 
00635         // write a tag to those sets denoting they're partition sets, with a value of the
00636         // proc number
00637         int* dum_ids = new int[nparts];  // this remains integer
00638         for( i = 0; i < nparts; i++ )
00639             dum_ids[i] = i;
00640 
00641         result = mbImpl->tag_set_data( part_set_tag, partSets, dum_ids );
00642         delete[] dum_ids;
00643 
00644         // assign entities to the relevant sets
00645         std::vector< EntityHandle > tmp_part_sets;
00646         std::copy( partSets.begin(), partSets.end(), std::back_inserter( tmp_part_sets ) );
00647         Range::iterator rit;
00648         for( i = 0, rit = elems.begin(); rit != elems.end(); rit++, i++ )
00649         {
00650             result = mbImpl->add_entities( tmp_part_sets[assignment[i]], &( *rit ), 1 );MB_CHK_ERR( result );
00651         }
00652 
00653         // check for empty sets, warn if there are any
00654         Range empty_sets;
00655         for( rit = partSets.begin(); rit != partSets.end(); rit++ )
00656         {
00657             int num_ents = 0;
00658             result       = mbImpl->get_number_entities_by_handle( *rit, num_ents );
00659             if( MB_SUCCESS != result || !num_ents ) empty_sets.insert( *rit );
00660         }
00661         if( !empty_sets.empty() )
00662         {
00663             std::cout << "WARNING: " << empty_sets.size() << " empty sets in partition: ";
00664             for( rit = empty_sets.begin(); rit != empty_sets.end(); rit++ )
00665                 std::cout << *rit << " ";
00666             std::cout << std::endl;
00667         }
00668     }
00669 
00670     if( write_as_tags )
00671     {
00672         if( sizeof( int ) != sizeof( idx_t ) )
00673         {
00674             // allocate idx_teger-size partitions
00675             // first we have to copy to int, then assign
00676             int* assg_int = new int[elems.size()];
00677             for( int k = 0; k < (int)elems.size(); k++ )
00678                 assg_int[k] = assignment[k];
00679             result = mbImpl->tag_set_data( part_set_tag, elems, assg_int );MB_CHK_ERR( result );
00680             delete[] assg_int;
00681         }
00682         else
00683             result = mbImpl->tag_set_data( part_set_tag, elems, assignment );MB_CHK_ERR( result );
00684     }
00685 
00686     return MB_SUCCESS;
00687 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines