MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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, 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" ) != 0 && strcmp( method, "ML_KWAY" ) != 0 ) 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 }