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