Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
ZoltanPartitioner.cpp
Go to the documentation of this file.
1 /* $Id$
2  *
3  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
4  * storing and accessing finite element mesh data.
5  *
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  *
10  * This library is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 2.1 of the License, or (at your option) any later version.
14  *
15  */
16 
17 /** Get a mesh from MOAB and write a Zoltan partition set back into MOAB and to
18  * a file
19  *
20  *
21  */
22 
23 #include <iostream>
24 #include <cassert>
25 #include <sstream>
26 #include <algorithm>
27 
28 #include "moab/ZoltanPartitioner.hpp"
29 #include "moab/Interface.hpp"
30 #include "Internals.hpp"
31 #include "moab/Range.hpp"
32 #include "moab/WriteUtilIface.hpp"
33 #include "moab/MeshTopoUtil.hpp"
34 #include "MBTagConventions.hpp"
35 #include "moab/CN.hpp"
36 // used for gnomonic projection
37 #include "moab/IntxMesh/IntxUtils.hpp"
38 
39 #ifdef MOAB_HAVE_CGM
40 #include "CGMConfig.h"
41 #include <limits>
42 #include "RefEntity.hpp"
43 #include "RefFace.hpp"
44 #include "RefEdge.hpp"
45 #include "RefVertex.hpp"
46 #include "Body.hpp"
47 #include "TopologyEntity.hpp"
48 #include "CastTo.hpp"
49 #include "CABodies.hpp"
50 #include "TDParallel.hpp"
51 #include "TDUniqueId.hpp"
52 #endif
53 
54 using namespace moab;
55 
56 #define RR \
57  if( MB_SUCCESS != result ) return result
58 
59 static double* Points = NULL;
60 static int* GlobalIds = NULL;
61 static int NumPoints = 0;
62 static int* NumEdges = NULL;
63 static int* NborGlobalId = NULL;
64 static int* NborProcs = NULL;
65 static double* ObjWeights = NULL;
66 static double* EdgeWeights = NULL;
67 static int* Parts = NULL;
68 
69 const bool debug = false;
70 
71 ZoltanPartitioner::ZoltanPartitioner( Interface* impl,
72 #ifdef MOAB_HAVE_MPI
73  ParallelComm* parcomm,
74 #endif
75  const bool use_coords,
76  int argc,
77  char** argv
78 #ifdef MOAB_HAVE_CGM
79  ,
80  GeometryQueryTool* gqt
81 #endif
82 
83  )
84  : PartitionerBase< int >( impl,
85  use_coords
86 #ifdef MOAB_HAVE_MPI
87  ,
88  parcomm
89 #endif
90  ),
91  myZZ( NULL ), myNumPts( 0 ), argcArg( argc ), argvArg( argv )
92 #ifdef MOAB_HAVE_CGM
93  ,
94  gti( gqt )
95 #endif
96 {
97 }
98 
99 ZoltanPartitioner::~ZoltanPartitioner()
100 {
101  if( NULL != myZZ ) delete myZZ;
102 }
103 
104 ErrorCode ZoltanPartitioner::balance_mesh( const char* zmethod,
105  const char* other_method,
106  const bool write_as_sets,
107  const bool write_as_tags )
108 {
109  if( !strcmp( zmethod, "RR" ) && !strcmp( zmethod, "RCB" ) && !strcmp( zmethod, "RIB" ) &&
110  !strcmp( zmethod, "HSFC" ) && !strcmp( zmethod, "Hypergraph" ) && !strcmp( zmethod, "PHG" ) &&
111  !strcmp( zmethod, "PARMETIS" ) && !strcmp( zmethod, "OCTPART" ) )
112  {
113  std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be "
114  << "RR, RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl;
115  return MB_FAILURE;
116  }
117 
118  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
119  std::vector< int > ids; // point ids from MOAB
120  std::vector< int > adjs, length;
121  Range elems;
122 
123  // Get a mesh from MOAB and divide it across processors.
124 
125  ErrorCode result;
126 
127  if( mbpc->proc_config().proc_rank() == 0 )
128  {
129  result = assemble_graph( 3, pts, ids, adjs, length, elems );RR;
130  }
131 
132  myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0] );
133 
134  // Initialize Zoltan. This is a C call. The simple C++ code
135  // that creates Zoltan objects does not keep track of whether
136  // Zoltan_Initialize has been called.
137 
138  float version;
139 
140  Zoltan_Initialize( argcArg, argvArg, &version );
141 
142  // Create Zoltan object. This calls Zoltan_Create.
143  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
144 
145  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
146  SetRCB_Parameters();
147  else if( !strcmp( zmethod, "RIB" ) )
148  SetRIB_Parameters();
149  else if( !strcmp( zmethod, "HSFC" ) )
150  SetHSFC_Parameters();
151  else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) )
152  if( NULL == other_method )
153  SetHypergraph_Parameters( "auto" );
154  else
155  SetHypergraph_Parameters( other_method );
156  else if( !strcmp( zmethod, "PARMETIS" ) )
157  {
158  if( NULL == other_method )
159  SetPARMETIS_Parameters( "RepartGDiffusion" );
160  else
161  SetPARMETIS_Parameters( other_method );
162  }
163  else if( !strcmp( zmethod, "OCTPART" ) )
164  {
165  if( NULL == other_method )
166  SetOCTPART_Parameters( "2" );
167  else
168  SetOCTPART_Parameters( other_method );
169  }
170 
171  // Call backs:
172 
173  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
174  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
175  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
176  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
177  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
178  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
179 
180  // Perform the load balancing partitioning
181 
182  int changes;
183  int numGidEntries;
184  int numLidEntries;
185  int numImport;
186  ZOLTAN_ID_PTR importGlobalIds;
187  ZOLTAN_ID_PTR importLocalIds;
188  int* importProcs;
189  int* importToPart;
190  int numExport;
191  ZOLTAN_ID_PTR exportGlobalIds;
192  ZOLTAN_ID_PTR exportLocalIds;
193  int* exportProcs;
194  int* exportToPart;
195 
196  int rc = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, numImport, importGlobalIds, importLocalIds,
197  importProcs, importToPart, numExport, exportGlobalIds, exportLocalIds, exportProcs,
198  exportToPart );
199 
200  rc = mbGlobalSuccess( rc );
201 
202  if( !rc )
203  {
204  mbPrintGlobalResult( "==============Result==============", myNumPts, numImport, numExport, changes );
205  }
206  else
207  {
208  return MB_FAILURE;
209  }
210 
211  // take results & write onto MOAB partition sets
212 
213  int* assignment;
214 
215  mbFinalizePoints( (int)ids.size(), numExport, exportLocalIds, exportProcs, &assignment );
216 
217  if( mbpc->proc_config().proc_rank() == 0 )
218  {
219  result = write_partition( mbpc->proc_config().proc_size(), elems, assignment, write_as_sets, write_as_tags );
220 
221  if( MB_SUCCESS != result ) return result;
222 
223  free( (int*)assignment );
224  }
225 
226  // Free the memory allocated for lists returned by LB_Parition()
227 
228  myZZ->LB_Free_Part( &importGlobalIds, &importLocalIds, &importProcs, &importToPart );
229  myZZ->LB_Free_Part( &exportGlobalIds, &exportLocalIds, &exportProcs, &exportToPart );
230 
231  // Implementation note: A Zoltan object contains an MPI communicator.
232  // When the Zoltan object is destroyed, it uses it's MPI communicator.
233  // So it is important that the Zoltan object is destroyed before
234  // the MPI communicator is destroyed. To ensure this, dynamically
235  // allocate the Zoltan object, so you can explicitly destroy it.
236  // If you create a Zoltan object on the stack, it's destructor will
237  // be invoked atexit, possibly after the communicator's
238  // destructor.
239 
240  return MB_SUCCESS;
241 }
242 
243 ErrorCode ZoltanPartitioner::repartition( std::vector< double >& x,
244  std::vector< double >& y,
245  std::vector< double >& z,
246  int StartID,
247  const char* zmethod,
248  Range& localGIDs )
249 {
250  //
251  int nprocs = mbpc->proc_config().proc_size();
252  int rank = mbpc->proc_config().proc_rank();
253  clock_t t = clock();
254  // form pts and ids, as in assemble_graph
255  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
256  pts.resize( x.size() * 3 );
257  std::vector< int > ids; // point ids from MOAB
258  ids.resize( x.size() );
259  for( size_t i = 0; i < x.size(); i++ )
260  {
261  pts[3 * i] = x[i];
262  pts[3 * i + 1] = y[i];
263  pts[3 * i + 2] = z[i];
264  ids[i] = StartID + (int)i;
265  }
266  // do not get out until done!
267 
268  Points = &pts[0];
269  GlobalIds = &ids[0];
270  NumPoints = (int)x.size();
271  NumEdges = NULL;
272  NborGlobalId = NULL;
273  NborProcs = NULL;
274  ObjWeights = NULL;
275  EdgeWeights = NULL;
276  Parts = NULL;
277 
278  float version;
279  if( rank == 0 ) std::cout << "Initializing zoltan..." << std::endl;
280 
281  Zoltan_Initialize( argcArg, argvArg, &version );
282 
283  // Create Zoltan object. This calls Zoltan_Create.
284  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
285 
286  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
287  SetRCB_Parameters();
288  else if( !strcmp( zmethod, "RIB" ) )
289  SetRIB_Parameters();
290  else if( !strcmp( zmethod, "HSFC" ) )
291  SetHSFC_Parameters();
292 
293  // set # requested partitions
294  char buff[10];
295  sprintf( buff, "%d", nprocs );
296  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
297  if( ZOLTAN_OK != retval ) return MB_FAILURE;
298 
299  // request all, import and export
300  retval = myZZ->Set_Param( "RETURN_LISTS", "ALL" );
301  if( ZOLTAN_OK != retval ) return MB_FAILURE;
302 
303  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
304  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
305  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
306  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
307  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
308  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
309 
310  // Perform the load balancing partitioning
311 
312  int changes;
313  int numGidEntries;
314  int numLidEntries;
315  int num_import;
316  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
317  int* import_procs;
318  int* import_to_part;
319  int num_export;
320  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
321  int *assign_procs, *assign_parts;
322 
323  if( rank == 0 )
324  std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nprocs
325  << " processors..." << std::endl;
326 
327  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids,
328  import_procs, import_to_part, num_export, export_global_ids, export_local_ids,
329  assign_procs, assign_parts );
330  if( ZOLTAN_OK != retval ) return MB_FAILURE;
331 
332  if( rank == 0 )
333  {
334  std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
335  t = clock();
336  }
337 
338  std::sort( import_global_ids, import_global_ids + num_import, std::greater< int >() );
339  std::sort( export_global_ids, export_global_ids + num_export, std::greater< int >() );
340 
341  Range iniGids( (EntityHandle)StartID, (EntityHandle)StartID + x.size() - 1 );
342  Range imported, exported;
343  std::copy( import_global_ids, import_global_ids + num_import, range_inserter( imported ) );
344  std::copy( export_global_ids, export_global_ids + num_export, range_inserter( exported ) );
345  localGIDs = subtract( iniGids, exported );
346  localGIDs = unite( localGIDs, imported );
347  // Free data structures allocated by Zoltan::LB_Partition
348  retval = myZZ->LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part );
349  if( ZOLTAN_OK != retval ) return MB_FAILURE;
350  retval = myZZ->LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts );
351  if( ZOLTAN_OK != retval ) return MB_FAILURE;
352 
353  return MB_SUCCESS;
354 }
355 
356 ErrorCode ZoltanPartitioner::partition_inferred_mesh( EntityHandle sfileset,
357  size_t num_parts,
358  int part_dim,
359  const bool write_as_sets,
360  int projection_type )
361 {
362  ErrorCode result;
363 
364  moab::Range elverts;
365  result = mbImpl->get_entities_by_dimension( sfileset, part_dim, elverts );RR;
366 
367  std::vector< double > elcoords( elverts.size() * 3 );
368  result = mbImpl->get_coords( elverts, &elcoords[0] );RR;
369 
370  std::vector< std::vector< EntityHandle > > part_assignments( num_parts );
371  int part, proc;
372 
373  // Loop over coordinates
374  for( size_t iel = 0; iel < elverts.size(); iel++ )
375  {
376  // Gather coordinates into temporary array
377  double* ecoords = &elcoords[iel * 3];
378 
379  // do a projection if needed
380  if( projection_type > 0 ) IntxUtils::transform_coordinates( ecoords, projection_type );
381 
382  // Compute the coordinate's part assignment
383  myZZ->LB_Point_PP_Assign( ecoords, proc, part );
384 
385  // Store the part assignment in the return array
386  part_assignments[part].push_back( elverts[iel] );
387  }
388 
389  // get the partition set tag
390  Tag part_set_tag = mbpc->partition_tag();
391 
392  // If some entity sets exist, let us delete those first.
393  if( write_as_sets )
394  {
395  moab::Range oldpsets;
396  result = mbImpl->get_entities_by_type_and_tag( sfileset, MBENTITYSET, &part_set_tag, NULL, 1, oldpsets,
397  Interface::UNION );RR;
398  result = mbImpl->remove_entities( sfileset, oldpsets );RR;
399  }
400 
401  size_t nparts_assigned = 0;
402  for( size_t partnum = 0; partnum < num_parts; ++partnum )
403  {
404  std::vector< EntityHandle >& partvec = part_assignments[partnum];
405 
406  nparts_assigned += ( partvec.size() ? 1 : 0 );
407 
408  if( write_as_sets )
409  {
410  EntityHandle partNset;
411  result = mbImpl->create_meshset( moab::MESHSET_SET, partNset );RR;
412  if( partvec.size() )
413  {
414  result = mbImpl->add_entities( partNset, &partvec[0], partvec.size() );RR;
415  }
416  result = mbImpl->add_parent_child( sfileset, partNset );RR;
417 
418  // part numbering should start from 0
419  int ipartnum = (int)partnum;
420 
421  // assign partitions as a sparse tag by grouping elements under sets
422  result = mbImpl->tag_set_data( part_set_tag, &partNset, 1, &ipartnum );RR;
423  }
424  else
425  {
426  /* assign as a dense vector to all elements
427  allocate integer-size partitions */
428  std::vector< int > assignment( partvec.size(),
429  partnum ); // initialize all values to partnum
430  result = mbImpl->tag_set_data( part_set_tag, &partvec[0], partvec.size(), &assignment[0] );RR;
431  }
432  }
433 
434  if( nparts_assigned != num_parts )
435  {
436  std::cout << "WARNING: The inference yielded lesser number of parts (" << nparts_assigned
437  << ") than requested by user (" << num_parts << ").\n";
438  }
439 
440  return MB_SUCCESS;
441 }
442 
443 ErrorCode ZoltanPartitioner::partition_mesh_and_geometry( const double part_geom_mesh_size,
444  const int nparts,
445  const char* zmethod,
446  const char* other_method,
447  double imbal_tol,
448  const int part_dim,
449  const bool write_as_sets,
450  const bool write_as_tags,
451  const int obj_weight,
452  const int edge_weight,
453 #ifdef MOAB_HAVE_CGM
454  const bool part_surf,
455  const bool ghost,
456 #else
457  const bool,
458  const bool,
459 #endif
460  const int projection_type,
461  const bool recompute_rcb_box,
462  const bool print_time )
463 {
464  // should only be called in serial
465  if( mbpc->proc_config().proc_size() != 1 )
466  {
467  std::cout << "ZoltanPartitioner::partition_mesh_and_geometry must be called in serial." << std::endl;
468  return MB_FAILURE;
469  }
470  clock_t t = clock();
471  if( NULL != zmethod && strcmp( zmethod, "RR" ) && strcmp( zmethod, "RCB" ) && strcmp( zmethod, "RIB" ) &&
472  strcmp( zmethod, "HSFC" ) && strcmp( zmethod, "Hypergraph" ) && strcmp( zmethod, "PHG" ) &&
473  strcmp( zmethod, "PARMETIS" ) && strcmp( zmethod, "OCTPART" ) )
474  {
475  std::cout << "ERROR node " << mbpc->proc_config().proc_rank() << ": Method must be "
476  << "RCB, RIB, HSFC, Hypergraph (PHG), PARMETIS, or OCTPART" << std::endl;
477  return MB_FAILURE;
478  }
479 
480  bool part_geom = false;
481  if( 0 == strcmp( zmethod, "RR" ) || 0 == strcmp( zmethod, "RCB" ) || 0 == strcmp( zmethod, "RIB" ) ||
482  0 == strcmp( zmethod, "HSFC" ) )
483  part_geom = true; // so no adjacency / edges needed
484  std::vector< double > pts; // x[0], y[0], z[0], ... from MOAB
485  std::vector< int > ids; // point ids from MOAB
486  std::vector< int > adjs, length, parts;
487  std::vector< double > obj_weights, edge_weights;
488  Range elems;
489 #ifdef MOAB_HAVE_CGM
490  DLIList< RefEntity* > entities;
491 #endif
492 
493  // Get a mesh from MOAB and diide it across processors.
494 
495  ErrorCode result = MB_SUCCESS;
496 
497  // short-circuit everything if RR partition is requested
498  if( !strcmp( zmethod, "RR" ) )
499  {
500  if( part_geom_mesh_size < 0. )
501  {
502  // get all elements
503  result = mbImpl->get_entities_by_dimension( 0, part_dim, elems );RR;
504  if( elems.empty() ) return MB_FAILURE;
505  // make a trivial assignment vector
506  std::vector< int > assign_vec( elems.size() );
507  int num_per = elems.size() / nparts;
508  int extra = elems.size() % nparts;
509  if( extra ) num_per++;
510  int nstart = 0;
511  for( int i = 0; i < nparts; i++ )
512  {
513  if( i == extra ) num_per--;
514  std::fill( &assign_vec[nstart], &assign_vec[nstart + num_per], i );
515  nstart += num_per;
516  }
517 
518  result = write_partition( nparts, elems, &assign_vec[0], write_as_sets, write_as_tags );RR;
519  }
520  else
521  {
522 #ifdef MOAB_HAVE_CGM
523  result = partition_round_robin( nparts );RR;
524 #endif
525  }
526 
527  return result;
528  }
529 
530  std::cout << "Assembling graph..." << std::endl;
531 
532  if( part_geom_mesh_size < 0. )
533  {
534  // if (!part_geom) {
535  result = assemble_graph( part_dim, pts, ids, adjs, length, elems, part_geom, projection_type );RR;
536  }
537  else
538  {
539 #ifdef MOAB_HAVE_CGM
540  result = assemble_graph( part_dim, pts, ids, adjs, length, obj_weights, edge_weights, parts, entities,
541  part_geom_mesh_size, nparts );RR;
542 
543  if( debug )
544  {
545  int n_ids = ids.size();
546  entities.reset();
547  int i_leng = 0;
548  for( int i = 0; i < n_ids; i++ )
549  {
550  std::cout << "graph_input_ids[" << i << "]=" << ids[i] << ",obj_weights=" << obj_weights[i]
551  << ",entity_id=" << entities.get_and_step()->id() << ",part=" << parts[i] << std::endl;
552  for( int j = 0; j < length[i]; j++ )
553  {
554  std::cout << "adjs[" << j << "]=" << adjs[i_leng] << ",edge_weights=" << edge_weights[i_leng]
555  << std::endl;
556  i_leng++;
557  }
558  }
559  }
560 #endif
561  }
562  if( print_time )
563  {
564  std::cout << " time to assemble graph: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
565  t = clock();
566  }
567  double* o_wgt = NULL;
568  double* e_wgt = NULL;
569  if( obj_weights.size() > 0 ) o_wgt = &obj_weights[0];
570  if( edge_weights.size() > 0 ) e_wgt = &edge_weights[0];
571 
572  myNumPts = mbInitializePoints( (int)ids.size(), &pts[0], &ids[0], &adjs[0], &length[0], o_wgt, e_wgt, &parts[0],
573  part_geom );
574 
575  // Initialize Zoltan. This is a C call. The simple C++ code
576  // that creates Zoltan objects does not keep track of whether
577  // Zoltan_Initialize has been called.
578 
579  if( print_time )
580  {
581  std::cout << " time to initialize points: " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
582  t = clock();
583  }
584  float version;
585 
586  std::cout << "Initializing zoltan..." << std::endl;
587 
588  Zoltan_Initialize( argcArg, argvArg, &version );
589 
590  // Create Zoltan object. This calls Zoltan_Create.
591  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
592 
593  if( NULL == zmethod || !strcmp( zmethod, "RCB" ) )
594  SetRCB_Parameters( recompute_rcb_box );
595  else if( !strcmp( zmethod, "RIB" ) )
596  SetRIB_Parameters();
597  else if( !strcmp( zmethod, "HSFC" ) )
598  SetHSFC_Parameters();
599  else if( !strcmp( zmethod, "Hypergraph" ) || !strcmp( zmethod, "PHG" ) )
600  {
601  if( NULL == other_method || ( other_method[0] == '\0' ) )
602  SetHypergraph_Parameters( "auto" );
603  else
604  SetHypergraph_Parameters( other_method );
605 
606  if( imbal_tol )
607  {
608  std::ostringstream str;
609  str << imbal_tol;
610  myZZ->Set_Param( "IMBALANCE_TOL", str.str().c_str() ); // no debug messages
611  }
612  }
613  else if( !strcmp( zmethod, "PARMETIS" ) )
614  {
615  if( NULL == other_method )
616  SetPARMETIS_Parameters( "RepartGDiffusion" );
617  else
618  SetPARMETIS_Parameters( other_method );
619  }
620  else if( !strcmp( zmethod, "OCTPART" ) )
621  {
622  if( NULL == other_method )
623  SetOCTPART_Parameters( "2" );
624  else
625  SetOCTPART_Parameters( other_method );
626  }
627 
628  // set # requested partitions
629  char buff[10];
630  sprintf( buff, "%d", nparts );
631  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
632  if( ZOLTAN_OK != retval ) return MB_FAILURE;
633 
634  // request only partition assignments
635  retval = myZZ->Set_Param( "RETURN_LISTS", "PARTITION ASSIGNMENTS" );
636  if( ZOLTAN_OK != retval ) return MB_FAILURE;
637 
638  if( obj_weight > 0 )
639  {
640  std::ostringstream str;
641  str << obj_weight;
642  retval = myZZ->Set_Param( "OBJ_WEIGHT_DIM", str.str().c_str() );
643  if( ZOLTAN_OK != retval ) return MB_FAILURE;
644  }
645 
646  if( edge_weight > 0 )
647  {
648  std::ostringstream str;
649  str << edge_weight;
650  retval = myZZ->Set_Param( "EDGE_WEIGHT_DIM", str.str().c_str() );
651  if( ZOLTAN_OK != retval ) return MB_FAILURE;
652  }
653 
654  // Call backs:
655 
656  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
657  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
658  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
659  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
660  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
661  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
662  if( part_geom_mesh_size > 0. )
663  {
664  myZZ->Set_Part_Multi_Fn( mbGetPart, NULL );
665  }
666 
667  // Perform the load balancing partitioning
668 
669  int changes;
670  int numGidEntries;
671  int numLidEntries;
672  int dumnum1;
673  ZOLTAN_ID_PTR dum_local, dum_global;
674  int *dum1, *dum2;
675  int num_assign;
676  ZOLTAN_ID_PTR assign_gid, assign_lid;
677  int *assign_procs, *assign_parts;
678 
679  std::cout << "Computing partition using " << ( zmethod ? zmethod : "RCB" ) << " method for " << nparts
680  << " processors..." << std::endl;
681 #ifndef NDEBUG
682 #if 0
683  if (NULL == zmethod || !strcmp(zmethod, "RCB"))
684  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 1, 0, 0);
685 
686  if ( !strcmp(zmethod, "PHG"))
687  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)zmethod, 1, 0, 1, 1);
688 #endif
689 #endif
690  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, dumnum1, dum_global, dum_local, dum1, dum2,
691  num_assign, assign_gid, assign_lid, assign_procs, assign_parts );
692  if( ZOLTAN_OK != retval ) return MB_FAILURE;
693 
694  if( print_time )
695  {
696  std::cout << " time to LB_partition " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
697  t = clock();
698  }
699  // take results & write onto MOAB partition sets
700  std::cout << "Saving partition information to MOAB..." << std::endl;
701 
702  if( part_geom_mesh_size < 0. )
703  {
704  // if (!part_geom) {
705  result = write_partition( nparts, elems, assign_parts, write_as_sets, write_as_tags );
706  }
707  else
708  {
709 #ifdef MOAB_HAVE_CGM
710  result = write_partition( nparts, entities, assign_parts, obj_weights, part_surf, ghost );
711 #endif
712  }
713 
714  if( print_time )
715  {
716  std::cout << " time to write partition in memory " << ( clock() - t ) / (double)CLOCKS_PER_SEC << "s. \n";
717  t = clock();
718  }
719 
720  if( MB_SUCCESS != result ) return result;
721 
722  // Free the memory allocated for lists returned by LB_Parition()
723  myZZ->LB_Free_Part( &assign_gid, &assign_lid, &assign_procs, &assign_parts );
724 
725  return MB_SUCCESS;
726 }
727 
728 ErrorCode ZoltanPartitioner::include_closure()
729 {
730  ErrorCode result;
731  Range ents;
732  Range adjs;
733  std::cout << "Adding closure..." << std::endl;
734 
735  for( Range::iterator rit = partSets.begin(); rit != partSets.end(); ++rit )
736  {
737 
738  // get the top-dimensional entities in the part
739  result = mbImpl->get_entities_by_handle( *rit, ents, true );RR;
740 
741  if( ents.empty() ) continue;
742 
743  // get intermediate-dimensional adjs and add to set
744  for( int d = mbImpl->dimension_from_handle( *ents.begin() ) - 1; d >= 1; d-- )
745  {
746  adjs.clear();
747  result = mbImpl->get_adjacencies( ents, d, false, adjs, Interface::UNION );RR;
748  result = mbImpl->add_entities( *rit, adjs );RR;
749  }
750 
751  // now get vertices and add to set; only need to do for ents, not for adjs
752  adjs.clear();
753  result = mbImpl->get_adjacencies( ents, 0, false, adjs, Interface::UNION );RR;
754  result = mbImpl->add_entities( *rit, adjs );RR;
755 
756  ents.clear();
757  }
758 
759  // now go over non-part entity sets, looking for contained entities
760  Range sets, part_ents;
761  result = mbImpl->get_entities_by_type( 0, MBENTITYSET, sets );RR;
762  for( Range::iterator rit = sets.begin(); rit != sets.end(); ++rit )
763  {
764  // skip parts
765  if( partSets.find( *rit ) != partSets.end() ) continue;
766 
767  // get entities in this set, recursively
768  ents.clear();
769  result = mbImpl->get_entities_by_handle( *rit, ents, true );RR;
770 
771  // now check over all parts
772  for( Range::iterator rit2 = partSets.begin(); rit2 != partSets.end(); ++rit2 )
773  {
774  part_ents.clear();
775  result = mbImpl->get_entities_by_handle( *rit2, part_ents, false );RR;
776  Range int_range = intersect( ents, part_ents );
777  if( !int_range.empty() )
778  {
779  // non-empty intersection, add to part set
780  result = mbImpl->add_entities( *rit2, &( *rit ), 1 );RR;
781  }
782  }
783  }
784 
785  // finally, mark all the part sets as having the closure
786  Tag closure_tag;
787  result =
788  mbImpl->tag_get_handle( "INCLUDES_CLOSURE", 1, MB_TYPE_INTEGER, closure_tag, MB_TAG_SPARSE | MB_TAG_CREAT );RR;
789 
790  std::vector< int > closure_vals( partSets.size(), 1 );
791  result = mbImpl->tag_set_data( closure_tag, partSets, &closure_vals[0] );RR;
792 
793  return MB_SUCCESS;
794 }
795 
796 ErrorCode ZoltanPartitioner::assemble_graph( const int dimension,
797  std::vector< double >& coords,
798  std::vector< int >& moab_ids,
799  std::vector< int >& adjacencies,
800  std::vector< int >& length,
801  Range& elems,
802  bool part_geom,
803  int projection_type )
804 {
805  // assemble a graph with vertices equal to elements of specified dimension, edges
806  // signified by list of other elements to which an element is connected
807 
808  // get the elements of that dimension
809  ErrorCode result = mbImpl->get_entities_by_dimension( 0, dimension, elems );
810  if( MB_SUCCESS != result || elems.empty() ) return result;
811 
812  // assign global ids
813  if( assign_global_ids )
814  {
815  EntityHandle rootset = 0;
816  result = mbpc->assign_global_ids( rootset, dimension, 1, true, true );RR;
817  }
818 
819  // now assemble the graph, calling MeshTopoUtil to get bridge adjacencies through d-1
820  // dimensional neighbors
821  MeshTopoUtil mtu( mbImpl );
822  Range adjs;
823  // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited
824  // by MBCN
825  int neighbors[5 * MAX_SUB_ENTITIES];
826  double avg_position[3];
827  int moab_id;
828 
829  // get the global id tag handle
830  Tag gid = mbImpl->globalId_tag();
831 
832  for( Range::iterator rit = elems.begin(); rit != elems.end(); ++rit )
833  {
834 
835  if( !part_geom )
836  {
837  // get bridge adjacencies
838  adjs.clear();
839  result = mtu.get_bridge_adjacencies( *rit, ( dimension > 0 ? dimension - 1 : 3 ), dimension, adjs );RR;
840 
841  // get the graph vertex ids of those
842  if( !adjs.empty() )
843  {
844  assert( adjs.size() < 5 * MAX_SUB_ENTITIES );
845  result = mbImpl->tag_get_data( gid, adjs, neighbors );RR;
846  }
847 
848  // copy those into adjacencies vector
849  length.push_back( (int)adjs.size() );
850  std::copy( neighbors, neighbors + adjs.size(), std::back_inserter( adjacencies ) );
851  }
852 
853  // get average position of vertices
854  result = mtu.get_average_position( *rit, avg_position );RR;
855 
856  // get the graph vertex id for this element
857  result = mbImpl->tag_get_data( gid, &( *rit ), 1, &moab_id );RR;
858 
859  // copy those into coords vector
860  moab_ids.push_back( moab_id );
861  // transform coordinates to spherical coordinates, if requested
862  if( projection_type > 0 ) IntxUtils::transform_coordinates( avg_position, projection_type );
863 
864  std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
865  }
866 
867  if( debug )
868  {
869  std::cout << "Length vector: " << std::endl;
870  std::copy( length.begin(), length.end(), std::ostream_iterator< int >( std::cout, ", " ) );
871  std::cout << std::endl;
872  std::cout << "Adjacencies vector: " << std::endl;
873  std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( std::cout, ", " ) );
874  std::cout << std::endl;
875  std::cout << "Moab_ids vector: " << std::endl;
876  std::copy( moab_ids.begin(), moab_ids.end(), std::ostream_iterator< int >( std::cout, ", " ) );
877  std::cout << std::endl;
878  std::cout << "Coords vector: " << std::endl;
879  std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( std::cout, ", " ) );
880  std::cout << std::endl;
881  }
882 
883  return MB_SUCCESS;
884 }
885 
886 #ifdef MOAB_HAVE_CGM
887 ErrorCode ZoltanPartitioner::assemble_graph( const int /* dimension */,
888  std::vector< double >& /* coords */,
889  std::vector< int >& moab_ids,
890  std::vector< int >& adjacencies,
891  std::vector< int >& length,
892  std::vector< double >& obj_weights,
893  std::vector< double >& edge_weights,
894  std::vector< int >& parts,
895  DLIList< RefEntity* >& entities,
896  const double part_geom_mesh_size,
897  const int n_part )
898 {
899  // get body vertex weights
900  DLIList< RefEntity* > body_list;
901  gti->ref_entity_list( "body", body_list, CUBIT_FALSE );
902  DLIList< RefFace* > all_shared_surfs;
903  int n_bodies = body_list.size();
904  int body_map_index = 1;
905  int surf_map_index = n_bodies + 1;
906  int n_bodies_proc = n_bodies / n_part; // # of entities per processor
907  int i_bodies_proc = n_bodies_proc; // entity index limit for each processor
908  int proc = 0;
909 
910  body_list.reset();
911  for( int i = 0; i < n_bodies; i++ )
912  {
913  // assign part in round-robin
914  if( i == i_bodies_proc )
915  {
916  proc++;
917  if( proc < n_part )
918  i_bodies_proc += n_bodies_proc;
919  else
920  {
921  proc %= n_part;
922  i_bodies_proc++;
923  }
924  }
925  parts.push_back( proc );
926 
927  // volume mesh generation load estimation
928  RefEntity* body = body_list.get_and_step();
929  double vertex_w = body->measure();
930  vertex_w /= part_geom_mesh_size * part_geom_mesh_size * part_geom_mesh_size;
931  vertex_w = 3.8559e-5 * vertex_w * log( vertex_w );
932 
933  // add child non-interface face weights
934  DLIList< RefFace* > shared_surfs;
935  DLIList< RefFace* > faces;
936  ( dynamic_cast< TopologyEntity* >( body ) )->ref_faces( faces );
937  int n_face = faces.size();
938  faces.reset();
939  for( int j = 0; j < n_face; j++ )
940  {
941  RefFace* face = faces.get_and_step();
942  TopologyEntity* te = CAST_TO( face, TopologyEntity );
943  if( te->bridge_manager()->number_of_bridges() > 1 )
944  { // shared
945  shared_surfs.append( face );
946  }
947  else
948  { // non-shared
949  vertex_w += estimate_face_mesh_load( face, face->measure() );
950  }
951  }
952 
953  int temp_index = body_map_index++;
954  body_vertex_map[body->id()] = temp_index;
955  moab_ids.push_back( temp_index ); // add body as graph vertex
956  obj_weights.push_back( vertex_w ); // add weight for body graph vertex
957 
958  if( debug )
959  {
960  std::cout << "body=" << body->id() << ",graph_id=" << temp_index << ",weight=" << vertex_w
961  << ",part=" << proc << std::endl;
962  }
963 
964  // treat shared surfaces connected to this body
965  int n_shared = shared_surfs.size();
966  shared_surfs.reset();
967  for( int k = 0; k < n_shared; k++ )
968  { // add adjacencies
969  RefFace* face = shared_surfs.get_and_step();
970  std::map< int, int >::iterator iter = surf_vertex_map.find( face->id() );
971  if( iter != surf_vertex_map.end() )
972  {
973  temp_index = ( *iter ).second;
974  }
975  else
976  {
977  temp_index = surf_map_index++;
978  surf_vertex_map[face->id()] = temp_index;
979  all_shared_surfs.append( face );
980  }
981  adjacencies.push_back( temp_index );
982  double tmp_sw = estimate_face_comm_load( face, part_geom_mesh_size );
983  edge_weights.push_back( tmp_sw );
984 
985  if( debug )
986  {
987  std::cout << "adjac=" << temp_index << ",weight=" << tmp_sw << std::endl;
988  }
989  }
990  length.push_back( n_shared );
991  }
992  entities += body_list;
993 
994  // add shared surface as graph vertex
995  int n_all_shared_surf = all_shared_surfs.size();
996  all_shared_surfs.reset();
997  for( int i = 0; i < n_all_shared_surf; i++ )
998  {
999  RefEntity* face = all_shared_surfs.get_and_step();
1000  moab_ids.push_back( surf_vertex_map[face->id()] );
1001  double face_mesh_load = estimate_face_mesh_load( face, part_geom_mesh_size );
1002  obj_weights.push_back( face_mesh_load );
1003 
1004  // set surface object graph edges
1005  int min_ind = -1;
1006  double min_load = std::numeric_limits< double >::max();
1007  ;
1008  DLIList< Body* > parents;
1009  dynamic_cast< TopologyEntity* >( face )->bodies( parents );
1010  int n_parents = parents.size();
1011  parents.reset();
1012  for( int j = 0; j < n_parents; j++ )
1013  {
1014  int body_gid = body_vertex_map[parents.get_and_step()->id()];
1015  adjacencies.push_back( body_gid ); // add adjacency
1016  edge_weights.push_back( estimate_face_comm_load( face, part_geom_mesh_size ) );
1017 
1018  if( obj_weights[body_gid - 1] < min_load )
1019  {
1020  min_ind = body_gid - 1;
1021  min_load = obj_weights[min_ind];
1022  }
1023  }
1024  length.push_back( n_parents );
1025  entities.append( dynamic_cast< RefEntity* >( face ) );
1026  obj_weights[min_ind] += face_mesh_load;
1027  parts.push_back( parts[min_ind] );
1028 
1029  if( debug )
1030  {
1031  std::cout << "shared_surf=" << face->id() << ",graph_id=" << surf_vertex_map[face->id()]
1032  << ",weight=" << face_mesh_load << ",part=" << parts[min_ind] << std::endl;
1033  }
1034  }
1035 
1036  for( size_t i = 0; i < obj_weights.size(); i++ )
1037  if( obj_weights[i] < 1. ) obj_weights[i] = 1.;
1038  for( size_t i = 0; i < edge_weights.size(); i++ )
1039  if( edge_weights[i] < 1. ) edge_weights[i] = 1.;
1040 
1041  return MB_SUCCESS;
1042 }
1043 
1044 double ZoltanPartitioner::estimate_face_mesh_load( RefEntity* face, const double h )
1045 {
1046  GeometryType type = CAST_TO( face, RefFace )->geometry_type();
1047  double n = face->measure() / h / h;
1048  double n_logn = n * log( n );
1049 
1050  if( type == PLANE_SURFACE_TYPE )
1051  {
1052  return 1.536168737505151e-4 * n_logn;
1053  }
1054  else if( type == SPLINE_SURFACE_TYPE )
1055  {
1056  return 5.910511018383144e-4 * n_logn;
1057  }
1058  else if( type == CONE_SURFACE_TYPE || type == SPHERE_SURFACE_TYPE || type == TORUS_SURFACE_TYPE )
1059  {
1060  return 2.352511671418708e-4 * n_logn;
1061  }
1062 
1063  return 0.0;
1064 }
1065 
1066 double ZoltanPartitioner::estimate_face_comm_load( RefEntity* face, const double h )
1067 {
1068  double peri = 0.;
1069 #if( ( CGM_MAJOR_VERSION == 14 && CGM_MINOR_VERSION > 2 ) || CGM_MAJOR_VERSION >= 15 )
1070  DLIList< DLIList< RefEdge* > > ref_edge_loops;
1071 #else
1072  DLIList< DLIList< RefEdge* >* > ref_edge_loops;
1073 #endif
1074  CAST_TO( face, RefFace )->ref_edge_loops( ref_edge_loops );
1075  ref_edge_loops.reset();
1076 
1077 #if( ( CGM_MAJOR_VERSION == 14 && CGM_MINOR_VERSION > 2 ) || CGM_MAJOR_VERSION >= 15 )
1078  for( int i = 0; i < ref_edge_loops.size(); i++ )
1079  {
1080  DLIList< RefEdge* > eloop = ref_edge_loops.get_and_step();
1081  eloop.reset();
1082  for( int j = 0; j < eloop.size(); j++ )
1083  {
1084  peri += eloop.get_and_step()->measure();
1085  }
1086  }
1087 #else
1088  for( int i = 0; i < ref_edge_loops.size(); i++ )
1089  {
1090  DLIList< RefEdge* >* eloop = ref_edge_loops.get_and_step();
1091  eloop->reset();
1092  for( int j = 0; j < eloop->size(); j++ )
1093  {
1094  peri += eloop->get_and_step()->measure();
1095  }
1096  }
1097 #endif
1098 
1099  // return 104*face->measure()/sqrt(3)/h/h + 56/3*peri/h;
1100  return ( 104 * face->measure() / sqrt( 3 ) / h / h + 56 / 3 * peri / h ) / 700000.;
1101 }
1102 
1103 ErrorCode ZoltanPartitioner::write_partition( const int nparts,
1104  DLIList< RefEntity* > entities,
1105  const int* assignment,
1106  std::vector< double >& obj_weights,
1107  const bool part_surf,
1108  const bool ghost )
1109 {
1110  ErrorCode result;
1111 
1112  // actuate CA_BODIES and turn on auto flag for other attributes
1113  CGMApp::instance()->attrib_manager()->register_attrib_type( CA_BODIES, "bodies", "BODIES", CABodies_creator,
1114  CUBIT_TRUE, CUBIT_TRUE, CUBIT_TRUE, CUBIT_TRUE,
1115  CUBIT_TRUE, CUBIT_FALSE );
1116  CGMApp::instance()->attrib_manager()->auto_flag( CUBIT_TRUE );
1117 
1118  // set partition info to Body at first
1119  int n_entities = entities.size();
1120  entities.reset();
1121  for( int i = 0; i < n_entities; i++ )
1122  {
1123  int proc = assignment[i];
1124  DLIList< int > shared_procs;
1125  RefEntity* entity = entities.get_and_step();
1126 
1127  if( entity->entity_type_info() == typeid( Body ) )
1128  {
1129  shared_procs.append( proc );
1130  TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1131  if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs );
1132 
1133  if( debug )
1134  {
1135  std::cout << "body" << entity->id() << "_is_partitioned_to_p" << proc << std::endl;
1136  }
1137 
1138  // assign to volumes, it should be removed in future
1139  DLIList< RefVolume* > volumes;
1140  ( dynamic_cast< TopologyEntity* >( entity ) )->ref_volumes( volumes );
1141  int n_vol = volumes.size();
1142  volumes.reset();
1143  for( int j = 0; j < n_vol; j++ )
1144  {
1145  RefEntity* vol = volumes.get_and_step();
1146  td_par = (TDParallel*)vol->get_TD( &TDParallel::is_parallel );
1147  if( td_par == NULL ) td_par = new TDParallel( vol, NULL, &shared_procs );
1148  }
1149  }
1150  }
1151 
1152  // set partition info to shared surfaces
1153  entities.reset();
1154  for( int i = 0; i < n_entities; i++ )
1155  {
1156  int proc = assignment[i];
1157  DLIList< int > shared_procs;
1158  RefEntity* entity = entities.get_and_step();
1159 
1160  if( entity->entity_type_info() == typeid( RefFace ) )
1161  { // surface
1162  DLIList< Body* > parents;
1163  ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parents );
1164  int n_parents = parents.size();
1165  if( n_parents != 2 )
1166  { // check # of parents
1167  std::cerr << "# of parent bodies of interface surface should be 2." << std::endl;
1168  return MB_FAILURE;
1169  }
1170  shared_procs.append( proc ); // local proc
1171  parents.reset();
1172  for( int j = 0; j < 2; j++ )
1173  {
1174  RefEntity* parent = parents.get_and_step();
1175  TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel );
1176 
1177  if( parent_td == NULL )
1178  {
1179  std::cerr << "parent body should have balanced process." << std::endl;
1180  return MB_FAILURE;
1181  }
1182  int temp_proc = parent_td->get_charge_proc();
1183  if( temp_proc != proc ) shared_procs.append( temp_proc ); // remote proc
1184  }
1185 
1186  if( shared_procs.size() > 1 )
1187  { // if it is shared by 2 processors
1188  int merge_id = TDUniqueId::get_unique_id( entity );
1189  TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1190  if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs, NULL, merge_id, 1 );
1191 
1192  if( debug )
1193  {
1194  std::cout << "surf" << entity->id() << "_is_partitioned_to_p";
1195  for( int j = 0; j < shared_procs.size(); j++ )
1196  {
1197  std::cout << "," << shared_procs[j];
1198  }
1199  std::cout << std::endl;
1200  }
1201  }
1202  }
1203  }
1204 
1205  // do non-shared surface partition too
1206  if( part_surf )
1207  {
1208  result = partition_surface( nparts, entities, assignment, obj_weights );RR;
1209  }
1210 
1211  // partition shared edges and vertex
1212  result = partition_child_entities( 1, nparts, part_surf, ghost );RR;
1213  result = partition_child_entities( 0, nparts, part_surf );RR;
1214 
1215  if( debug )
1216  {
1217  entities.reset();
1218  for( int i = 0; i < n_entities; i++ )
1219  {
1220  RefEntity* entity = entities.get_and_step();
1221  if( entity->entity_type_info() == typeid( Body ) )
1222  {
1223  TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1224  CubitString st = entity->entity_name();
1225  DLIList< int >* sp = td_par->get_shared_proc_list();
1226  int n_sp = sp->size();
1227  sp->reset();
1228  for( int j = 0; j < n_sp; j++ )
1229  {
1230  std::cout << "partitioned_" << st.c_str() << ",proc=" << sp->get_and_step() << std::endl;
1231  }
1232  DLIList< int >* gp = td_par->get_ghost_proc_list();
1233  int n_gp = gp->size();
1234  sp->reset();
1235  for( int j = 0; j < n_gp; j++ )
1236  {
1237  std::cout << "partitioned_" << st.c_str() << ",ghost=" << gp->get_and_step() << std::endl;
1238  }
1239  }
1240  }
1241  }
1242 
1243  return MB_SUCCESS;
1244 }
1245 
1246 ErrorCode ZoltanPartitioner::partition_surface( const int nparts,
1247  DLIList< RefEntity* > entities,
1248  const int* assignment,
1249  std::vector< double >& obj_weights )
1250 {
1251  int i;
1252  double ave_load = 0.;
1253  double* loads = new double[nparts];
1254  for( i = 0; i < nparts; i++ )
1255  loads[i] = 0.;
1256 
1257  int n_entities = entities.size();
1258  entities.reset();
1259  for( i = 0; i < n_entities; i++ )
1260  {
1261  loads[assignment[i]] += obj_weights[i];
1262  ave_load += obj_weights[i];
1263  }
1264  ave_load /= nparts;
1265 
1266  if( debug )
1267  {
1268  for( i = 0; i < nparts; i++ )
1269  {
1270  std::cout << "loads_before_surface_partition[" << i << "]=" << loads[i] << std::endl;
1271  }
1272  }
1273 
1274  int n_iter = 0;
1275  bool b_stop = false;
1276  do
1277  {
1278  // get max and min load processors
1279  int max_proc = nparts, min_proc = 0;
1280  double min_load = std::numeric_limits< double >::max();
1281  double max_load = std::numeric_limits< double >::min();
1282  for( i = 0; i < nparts; i++ )
1283  {
1284  if( loads[i] < min_load )
1285  {
1286  min_load = loads[i];
1287  min_proc = i;
1288  }
1289  if( loads[i] > max_load )
1290  {
1291  max_load = loads[i];
1292  max_proc = i;
1293  }
1294  }
1295 
1296  double load_diff = max_load - ave_load;
1297  if( load_diff > ave_load / 10. )
1298  {
1299  bool b_moved = false;
1300  entities.reset();
1301  for( i = 0; i < n_entities; i++ )
1302  {
1303  if( b_moved ) break;
1304  if( assignment[i] == max_proc && // find maximum load processor bodies
1305  entities[i]->entity_type_info() == typeid( Body ) )
1306  {
1307  DLIList< RefFace* > faces;
1308  ( dynamic_cast< TopologyEntity* >( entities[i] ) )->ref_faces( faces );
1309  int n_face = faces.size();
1310  faces.reset();
1311  for( int j = 0; j < n_face; j++ )
1312  {
1313  RefFace* face = faces.get_and_step();
1314  TopologyEntity* te = CAST_TO( face, TopologyEntity );
1315  if( te->bridge_manager()->number_of_bridges() < 2 )
1316  { // non-shared
1317  TDParallel* td_par = (TDParallel*)face->get_TD( &TDParallel::is_parallel );
1318  if( td_par == NULL )
1319  { // only consider unpartitioned surfaces
1320  double face_load = face->measure();
1321  if( load_diff > min_load + face_load - ave_load )
1322  {
1323  loads[max_proc] -= face_load;
1324  loads[min_proc] += face_load;
1325  int merge_id = TDUniqueId::get_unique_id( face );
1326  DLIList< int > shared_procs;
1327  shared_procs.append( min_proc );
1328  shared_procs.append( max_proc );
1329  td_par = new TDParallel( face, NULL, &shared_procs, NULL, merge_id, 1 );
1330 
1331  if( debug )
1332  {
1333  std::cout << "non-shared surface " << face->id() << " is moved from p"
1334  << max_proc << " to p" << min_proc << std::endl;
1335  }
1336  b_moved = true;
1337  break;
1338  }
1339  }
1340  }
1341  }
1342  }
1343  }
1344  }
1345  else
1346  b_stop = true;
1347 
1348  n_iter++;
1349  } while( !b_stop && n_iter < 50 );
1350 
1351  if( debug )
1352  {
1353  for( i = 0; i < nparts; i++ )
1354  {
1355  std::cout << "loads_after_surface_partition[" << i << "]=" << loads[i] << std::endl;
1356  }
1357  }
1358 
1359  delete loads;
1360  return MB_SUCCESS;
1361 }
1362 
1363 ErrorCode ZoltanPartitioner::partition_round_robin( const int n_part )
1364 {
1365  int i, j, k;
1366  double* loads = new double[n_part]; // estimated loads for each processor
1367  double* ve_loads = new double[n_part]; // estimated loads for each processor
1368  for( i = 0; i < n_part; i++ )
1369  {
1370  loads[i] = 0.0;
1371  ve_loads[i] = 0.0;
1372  }
1373  DLIList< RefEntity* > body_entity_list;
1374  gti->ref_entity_list( "body", body_entity_list, CUBIT_FALSE );
1375  int n_entity = body_entity_list.size();
1376  int n_entity_proc = n_entity / n_part; // # of entities per processor
1377  int i_entity_proc = n_entity_proc; // entity index limit for each processor
1378  int proc = 0;
1379  RefEntity* entity;
1380 
1381  // assign processors to bodies
1382  body_entity_list.reset();
1383  for( i = 0; i < n_entity; i++ )
1384  {
1385  if( i == i_entity_proc )
1386  {
1387  proc++;
1388  if( proc < n_part )
1389  i_entity_proc += n_entity_proc;
1390  else
1391  {
1392  proc %= n_part;
1393  i_entity_proc++;
1394  }
1395  }
1396 
1397  // assign to bodies
1398  entity = body_entity_list.get_and_step();
1399  DLIList< int > shared_procs;
1400  shared_procs.append( proc );
1401  TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1402  if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs );
1403  loads[proc] += entity->measure();
1404 
1405  // assign to volumes, it should be removed in future
1406  DLIList< RefVolume* > volumes;
1407  ( dynamic_cast< TopologyEntity* >( entity ) )->ref_volumes( volumes );
1408  int n_vol = volumes.size();
1409  volumes.reset();
1410  for( j = 0; j < n_vol; j++ )
1411  {
1412  RefEntity* vol = volumes.get_and_step();
1413  td_par = (TDParallel*)vol->get_TD( &TDParallel::is_parallel );
1414  if( td_par == NULL ) td_par = new TDParallel( vol, NULL, &shared_procs );
1415  }
1416 
1417  // add local surface load
1418  DLIList< RefFace* > faces;
1419  ( dynamic_cast< TopologyEntity* >( entity ) )->ref_faces( faces );
1420  int n_face = faces.size();
1421  faces.reset();
1422  for( j = 0; j < n_face; j++ )
1423  {
1424  RefFace* face = faces.get_and_step();
1425  TopologyEntity* te = CAST_TO( face, TopologyEntity );
1426  if( te->bridge_manager()->number_of_bridges() < 2 ) loads[proc] = loads[proc] + face->measure();
1427  }
1428 
1429  // Get all child entities
1430  DLIList< RefEntity* > child_list;
1431  RefEntity::get_all_child_ref_entities( body_entity_list, child_list );
1432  int n_child = child_list.size();
1433 
1434  // assign processors to interface entities
1435  child_list.reset();
1436  for( j = 0; j < n_child; j++ )
1437  {
1438  entity = child_list.get_and_step();
1439  TopologyEntity* te = CAST_TO( entity, TopologyEntity );
1440 
1441  if( te->bridge_manager()->number_of_bridges() > 1 )
1442  {
1443  DLIList< Body* > parent_bodies;
1444  DLIList< int > child_shared_procs; // Shared processors of each child entity
1445  ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parent_bodies );
1446  int n_parent = parent_bodies.size();
1447 
1448  for( k = 0; k < n_parent; k++ )
1449  {
1450  RefEntity* parent_vol = CAST_TO( parent_bodies.get_and_step(), RefEntity );
1451  TDParallel* parent_td = (TDParallel*)parent_vol->get_TD( &TDParallel::is_parallel );
1452 
1453  if( parent_td == NULL )
1454  {
1455  PRINT_ERROR( "parent Volume has to be partitioned." );
1456  delete[] loads;
1457  delete[] ve_loads;
1458  return MB_FAILURE;
1459  }
1460  child_shared_procs.append_unique( parent_td->get_charge_proc() );
1461  }
1462 
1463  if( child_shared_procs.size() > 1 )
1464  { // if it is interface
1465  td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1466  if( td_par == NULL )
1467  {
1468  int merge_id = TDUniqueId::get_unique_id( entity );
1469  if( entity->entity_type_info() == typeid( RefFace ) )
1470  { // face
1471  if( child_shared_procs.size() != 2 )
1472  {
1473  PRINT_ERROR( "Error: # of shared processors of interface surface "
1474  "should be 2." );
1475  delete[] loads;
1476  delete[] ve_loads;
1477  return MB_FAILURE;
1478  }
1479 
1480  // balance interface surface loads
1481  if( loads[child_shared_procs[0]] > loads[child_shared_procs[1]] )
1482  child_shared_procs.reverse();
1483 
1484  loads[child_shared_procs[0]] = loads[child_shared_procs[0]] + entity->measure();
1485  td_par = new TDParallel( entity, NULL, &child_shared_procs, NULL, merge_id, 1 );
1486  } // face
1487  else if( entity->entity_type_info() == typeid( RefEdge ) ||
1488  entity->entity_type_info() == typeid( RefVertex ) )
1489  { // edge or vertex
1490  // balance interface surface loads
1491  int min_p = child_shared_procs[0];
1492  int n_shared_proc = child_shared_procs.size();
1493  for( k = 1; k < n_shared_proc; k++ )
1494  {
1495  if( ve_loads[child_shared_procs[k]] < ve_loads[min_p] ) min_p = child_shared_procs[k];
1496  }
1497  ve_loads[min_p] = ve_loads[min_p] + entity->measure();
1498  child_shared_procs.remove( min_p );
1499  child_shared_procs.insert_first( min_p );
1500  td_par = new TDParallel( entity, NULL, &child_shared_procs, NULL, merge_id, 1 );
1501  } // edge or vertex
1502  } // if (td_par == NULL)
1503  } // if it is interface
1504  } // if (te->bridge_manager()->number_of_bridges() > 1)
1505  } // for (j = 0; j < n_child; j++)
1506  } // for (i = 0; i < n_entity; i++)
1507 
1508  delete[] loads;
1509  delete[] ve_loads;
1510 
1511  return MB_SUCCESS;
1512 }
1513 
1514 // partition child entities to one of parent entity shared processors
1515 ErrorCode ZoltanPartitioner::partition_child_entities( const int dim,
1516  const int n_part,
1517  const bool part_surf,
1518  const bool ghost )
1519 {
1520  DLIList< RefEntity* > entity_list;
1521  if( dim == 0 )
1522  gti->ref_entity_list( "vertex", entity_list, CUBIT_FALSE );
1523  else if( dim == 1 )
1524  gti->ref_entity_list( "curve", entity_list, CUBIT_FALSE );
1525  else
1526  {
1527  std::cerr << "Dimention should be from 0 to 1." << std::endl;
1528  return MB_FAILURE;
1529  }
1530 
1531  int i, j, k;
1532  int n_entity = entity_list.size();
1533  double* loads = new double[n_part];
1534  for( i = 0; i < n_part; i++ )
1535  loads[i] = 0.;
1536  entity_list.reset();
1537 
1538  for( i = 0; i < n_entity; i++ )
1539  { // for all entities
1540  RefEntity* entity = entity_list.get_and_step();
1541  TopologyEntity* te = CAST_TO( entity, TopologyEntity );
1542 
1543  if( !part_surf && te->bridge_manager()->number_of_bridges() < 2 ) continue;
1544 
1545  DLIList< int > shared_procs;
1546  DLIList< Body* > parents;
1547  ( dynamic_cast< TopologyEntity* >( entity ) )->bodies( parents );
1548  int n_parents = parents.size();
1549  std::set< int > s_proc;
1550  parents.reset();
1551 
1552  // get shared processors from parent bodies
1553  for( j = 0; j < n_parents; j++ )
1554  {
1555  RefEntity* parent = parents.get_and_step();
1556  TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel );
1557 
1558  if( parent_td != NULL )
1559  {
1560  DLIList< int >* parent_procs = parent_td->get_shared_proc_list();
1561  int n_shared = parent_procs->size();
1562  parent_procs->reset();
1563  for( k = 0; k < n_shared; k++ )
1564  {
1565  int p = parent_procs->get_and_step();
1566  s_proc.insert( p );
1567  }
1568  }
1569  }
1570 
1571  if( part_surf )
1572  { // also get shared procs from parent surfaces
1573  DLIList< RefFace* > parent_faces;
1574  ( dynamic_cast< TopologyEntity* >( entity ) )->ref_faces( parent_faces );
1575  int n_pface = parent_faces.size();
1576  parent_faces.reset();
1577 
1578  // get shared processors from parent faces
1579  for( j = 0; j < n_pface; j++ )
1580  {
1581  RefEntity* parent = parent_faces.get_and_step();
1582  TDParallel* parent_td = (TDParallel*)parent->get_TD( &TDParallel::is_parallel );
1583 
1584  if( parent_td != NULL )
1585  {
1586  DLIList< int >* parent_procs = parent_td->get_shared_proc_list();
1587  int n_shared = parent_procs->size();
1588  parent_procs->reset();
1589 
1590  for( k = 0; k < n_shared; k++ )
1591  {
1592  int p = parent_procs->get_and_step();
1593  s_proc.insert( p );
1594  }
1595  }
1596  }
1597  }
1598 
1599  // find the minimum load processor and put it
1600  // at the front of the shared_procs list
1601  if( s_proc.size() > 1 )
1602  {
1603  int min_proc = 0;
1604  double min_load = std::numeric_limits< double >::max();
1605  std::set< int >::iterator iter = s_proc.begin();
1606  std::set< int >::iterator end_iter = s_proc.end();
1607  for( ; iter != end_iter; ++iter )
1608  {
1609  if( loads[*iter] < min_load )
1610  {
1611  min_load = loads[*iter];
1612  min_proc = *iter;
1613  }
1614  }
1615 
1616  if( dim == 1 )
1617  loads[min_proc] += entity->measure();
1618  else if( dim == 0 )
1619  loads[min_proc] += 1.;
1620  shared_procs.append( min_proc );
1621  iter = s_proc.begin();
1622  end_iter = s_proc.end();
1623  for( ; iter != end_iter; ++iter )
1624  {
1625  if( *iter != min_proc )
1626  {
1627  shared_procs.append( *iter );
1628  }
1629  }
1630 
1631  // add ghost geometries to shared processors for edge
1632  if( ghost )
1633  {
1634  parents.reset();
1635  for( j = 0; j < n_parents; j++ )
1636  { // for all parent bodies
1637  RefEntity* parent_vol = CAST_TO( parents.get_and_step(), RefEntity );
1638  TDParallel* parent_td = (TDParallel*)parent_vol->get_TD( &TDParallel::is_parallel );
1639  int n_shared_proc = shared_procs.size();
1640 
1641  for( k = 0; k < n_shared_proc; k++ )
1642  {
1643  parent_td->add_ghost_proc( shared_procs[k] );
1644  }
1645  }
1646  }
1647 
1648  // initialize tool data
1649  int merge_id = TDUniqueId::get_unique_id( entity );
1650  TDParallel* td_par = (TDParallel*)entity->get_TD( &TDParallel::is_parallel );
1651  if( td_par == NULL ) td_par = new TDParallel( entity, NULL, &shared_procs, NULL, merge_id, 1 );
1652  }
1653  }
1654 
1655  delete loads;
1656  return MB_SUCCESS;
1657 }
1658 #endif
1659 
1660 ErrorCode ZoltanPartitioner::write_partition( const int nparts,
1661  Range& elems,
1662  const int* assignment,
1663  const bool write_as_sets,
1664  const bool write_as_tags )
1665 {
1666  ErrorCode result;
1667 
1668  // get the partition set tag
1669  Tag part_set_tag;
1670  int dum_id = -1, i;
1671  result = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_set_tag,
1672  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );RR;
1673 
1674  // get any sets already with this tag, and clear them
1675  Range tagged_sets;
1676  result =
1677  mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &part_set_tag, NULL, 1, tagged_sets, Interface::UNION );RR;
1678  if( !tagged_sets.empty() )
1679  {
1680  result = mbImpl->clear_meshset( tagged_sets );RR;
1681  if( !write_as_sets )
1682  {
1683  result = mbImpl->tag_delete_data( part_set_tag, tagged_sets );RR;
1684  }
1685  }
1686 
1687  if( write_as_sets )
1688  {
1689  // first, create partition sets and store in vector
1690  partSets.clear();
1691 
1692  if( nparts > (int)tagged_sets.size() )
1693  {
1694  // too few partition sets - create missing ones
1695  int num_new = nparts - tagged_sets.size();
1696  for( i = 0; i < num_new; i++ )
1697  {
1698  EntityHandle new_set;
1699  result = mbImpl->create_meshset( MESHSET_SET, new_set );RR;
1700  tagged_sets.insert( new_set );
1701  }
1702  }
1703  else if( nparts < (int)tagged_sets.size() )
1704  {
1705  // too many partition sets - delete extras
1706  int num_del = tagged_sets.size() - nparts;
1707  for( i = 0; i < num_del; i++ )
1708  {
1709  EntityHandle old_set = tagged_sets.pop_back();
1710  result = mbImpl->delete_entities( &old_set, 1 );RR;
1711  }
1712  }
1713  // assign partition sets to vector
1714  partSets.swap( tagged_sets );
1715 
1716  // write a tag to those sets denoting they're partition sets, with a value of the
1717  // proc number
1718  int* dum_ids = new int[nparts];
1719  for( i = 0; i < nparts; i++ )
1720  dum_ids[i] = i;
1721 
1722  result = mbImpl->tag_set_data( part_set_tag, partSets, dum_ids );RR;
1723  // found out by valgrind when we run mbpart
1724  delete[] dum_ids;
1725  dum_ids = NULL;
1726 
1727  // assign entities to the relevant sets
1728  std::vector< EntityHandle > tmp_part_sets;
1729  // int N = (int)elems.size();
1730  std::copy( partSets.begin(), partSets.end(), std::back_inserter( tmp_part_sets ) );
1731  /*Range::reverse_iterator riter;
1732  for (i = N-1, riter = elems.rbegin(); riter != elems.rend(); ++riter, i--) {
1733  int assigned_part = assignment[i];
1734  part_ranges[assigned_part].insert(*riter);
1735  //result = mbImpl->add_entities(tmp_part_sets[assignment[i]], &(*rit), 1); RR;
1736  }*/
1737 
1738  Range::iterator rit;
1739  for( i = 0, rit = elems.begin(); rit != elems.end(); ++rit, i++ )
1740  {
1741  result = mbImpl->add_entities( tmp_part_sets[assignment[i]], &( *rit ), 1 );RR;
1742  }
1743  /*for (i=0; i<nparts; i++)
1744  {
1745  result = mbImpl->add_entities(tmp_part_sets[i], part_ranges[i]); RR;
1746  }
1747  delete [] part_ranges;*/
1748  // check for empty sets, warn if there are any
1749  Range empty_sets;
1750  for( rit = partSets.begin(); rit != partSets.end(); ++rit )
1751  {
1752  int num_ents = 0;
1753  result = mbImpl->get_number_entities_by_handle( *rit, num_ents );
1754  if( MB_SUCCESS != result || !num_ents ) empty_sets.insert( *rit );
1755  }
1756  if( !empty_sets.empty() )
1757  {
1758  std::cout << "WARNING: " << empty_sets.size() << " empty sets in partition: ";
1759  for( rit = empty_sets.begin(); rit != empty_sets.end(); ++rit )
1760  std::cout << *rit << " ";
1761  std::cout << std::endl;
1762  }
1763  }
1764 
1765  if( write_as_tags )
1766  {
1767  // allocate integer-size partitions
1768  result = mbImpl->tag_set_data( part_set_tag, elems, assignment );RR;
1769  }
1770 
1771  return MB_SUCCESS;
1772 }
1773 
1774 void ZoltanPartitioner::SetRCB_Parameters( const bool recompute_rcb_box )
1775 {
1776  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nRecursive Coordinate Bisection" << std::endl;
1777  // General parameters:
1778 
1779  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1780  myZZ->Set_Param( "LB_METHOD", "RCB" ); // recursive coordinate bisection
1781  // myZZ->Set_Param( "RCB_RECOMPUTE_BOX", "1" ); // recompute RCB box if needed ?
1782 
1783  // RCB parameters:
1784 
1785  myZZ->Set_Param( "RCB_OUTPUT_LEVEL", "1" );
1786  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition so that we can infer partitions
1787  // myZZ->Set_Param("RCB_RECTILINEAR_BLOCKS", "1"); // don't split point on boundary
1788  if( recompute_rcb_box ) myZZ->Set_Param( "RCB_RECOMPUTE_BOX", "1" );
1789 }
1790 
1791 void ZoltanPartitioner::SetRIB_Parameters()
1792 {
1793  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nRecursive Inertial Bisection" << std::endl;
1794  // General parameters:
1795 
1796  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1797  myZZ->Set_Param( "LB_METHOD", "RIB" ); // Recursive Inertial Bisection
1798 
1799  // RIB parameters:
1800 
1801  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition
1802  myZZ->Set_Param( "AVERAGE_CUTS", "1" );
1803 }
1804 
1805 void ZoltanPartitioner::SetHSFC_Parameters()
1806 {
1807  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nHilbert Space Filling Curve" << std::endl;
1808  // General parameters:
1809 
1810  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1811  myZZ->Set_Param( "LB_METHOD", "HSFC" ); // perform Hilbert space filling curve
1812 
1813  // HSFC parameters:
1814 
1815  myZZ->Set_Param( "KEEP_CUTS", "1" ); // save decomposition
1816 }
1817 
1818 void ZoltanPartitioner::SetHypergraph_Parameters( const char* phg_method )
1819 {
1820  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nHypergraph (or PHG): " << std::endl;
1821  // General parameters:
1822 
1823  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1824  myZZ->Set_Param( "LB_METHOD", "Hypergraph" ); // Hypergraph (or PHG)
1825 
1826  // Hypergraph (or PHG) parameters:
1827  myZZ->Set_Param( "PHG_COARSEPARTITION_METHOD", phg_method ); // CoarsePartitionMethod
1828 }
1829 
1830 void ZoltanPartitioner::SetPARMETIS_Parameters( const char* parmetis_method )
1831 {
1832  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nPARMETIS: " << parmetis_method << std::endl;
1833  // General parameters:
1834 
1835  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1836  myZZ->Set_Param( "LB_METHOD", "PARMETIS" ); // the ParMETIS library
1837 
1838  // PARMETIS parameters:
1839 
1840  myZZ->Set_Param( "PARMETIS_METHOD", parmetis_method ); // method in the library
1841 }
1842 
1843 void ZoltanPartitioner::SetOCTPART_Parameters( const char* oct_method )
1844 {
1845  if( mbpc->proc_config().proc_rank() == 0 ) std::cout << "\nOctree Partitioning: " << oct_method << std::endl;
1846  // General parameters:
1847 
1848  myZZ->Set_Param( "DEBUG_LEVEL", "0" ); // no debug messages
1849  myZZ->Set_Param( "LB_METHOD", "OCTPART" ); // octree partitioning
1850 
1851  // OCTPART parameters:
1852 
1853  myZZ->Set_Param( "OCT_METHOD", oct_method ); // the SFC to be used
1854  myZZ->Set_Param( "OCT_OUTPUT_LEVEL", "3" );
1855 }
1856 
1857 int ZoltanPartitioner::mbInitializePoints( int npts,
1858  double* pts,
1859  int* ids,
1860  int* adjs,
1861  int* length,
1862  double* obj_weights,
1863  double* edge_weights,
1864  int* parts,
1865  bool part_geom )
1866 {
1867  unsigned int i;
1868  int j;
1869  int *numPts, *nborProcs = NULL;
1870  int sum, ptsPerProc, ptsAssigned, mySize;
1871  MPI_Status stat;
1872  double* sendPts;
1873  int* sendIds;
1874  int* sendEdges = NULL;
1875  int* sendNborId = NULL;
1876  int* sendProcs;
1877 
1878  if( mbpc->proc_config().proc_rank() == 0 )
1879  {
1880  /* divide pts to start */
1881 
1882  numPts = (int*)malloc( sizeof( int ) * mbpc->proc_config().proc_size() );
1883  ptsPerProc = npts / mbpc->proc_config().proc_size();
1884  ptsAssigned = 0;
1885 
1886  for( i = 0; i < mbpc->proc_config().proc_size() - 1; i++ )
1887  {
1888  numPts[i] = ptsPerProc;
1889  ptsAssigned += ptsPerProc;
1890  }
1891 
1892  numPts[mbpc->proc_config().proc_size() - 1] = npts - ptsAssigned;
1893 
1894  mySize = numPts[mbpc->proc_config().proc_rank()];
1895  sendPts = pts + ( 3 * numPts[0] );
1896  sendIds = ids + numPts[0];
1897  sum = 0; // possible no adjacency sent
1898  if( !part_geom )
1899  {
1900  sendEdges = length + numPts[0];
1901 
1902  for( j = 0; j < numPts[0]; j++ )
1903  sum += length[j];
1904 
1905  sendNborId = adjs + sum;
1906 
1907  for( j = numPts[0]; j < npts; j++ )
1908  sum += length[j];
1909 
1910  nborProcs = (int*)malloc( sizeof( int ) * sum );
1911  }
1912  for( j = 0; j < sum; j++ )
1913  if( ( i = adjs[j] / ptsPerProc ) < mbpc->proc_config().proc_size() )
1914  nborProcs[j] = i;
1915  else
1916  nborProcs[j] = mbpc->proc_config().proc_size() - 1;
1917 
1918  sendProcs = nborProcs + ( sendNborId - adjs );
1919 
1920  for( i = 1; i < mbpc->proc_config().proc_size(); i++ )
1921  {
1922  MPI_Send( &numPts[i], 1, MPI_INT, i, 0x00, mbpc->comm() );
1923  MPI_Send( sendPts, 3 * numPts[i], MPI_DOUBLE, i, 0x01, mbpc->comm() );
1924  MPI_Send( sendIds, numPts[i], MPI_INT, i, 0x03, mbpc->comm() );
1925  MPI_Send( sendEdges, numPts[i], MPI_INT, i, 0x06, mbpc->comm() );
1926  sum = 0;
1927 
1928  for( j = 0; j < numPts[i]; j++ )
1929  sum += sendEdges[j];
1930 
1931  MPI_Send( sendNborId, sum, MPI_INT, i, 0x07, mbpc->comm() );
1932  MPI_Send( sendProcs, sum, MPI_INT, i, 0x08, mbpc->comm() );
1933  sendPts += ( 3 * numPts[i] );
1934  sendIds += numPts[i];
1935  sendEdges += numPts[i];
1936  sendNborId += sum;
1937  sendProcs += sum;
1938  }
1939 
1940  free( numPts );
1941  }
1942  else
1943  {
1944  MPI_Recv( &mySize, 1, MPI_INT, 0, 0x00, mbpc->comm(), &stat );
1945  pts = (double*)malloc( sizeof( double ) * 3 * mySize );
1946  ids = (int*)malloc( sizeof( int ) * mySize );
1947  length = (int*)malloc( sizeof( int ) * mySize );
1948  if( obj_weights != NULL ) obj_weights = (double*)malloc( sizeof( double ) * mySize );
1949  MPI_Recv( pts, 3 * mySize, MPI_DOUBLE, 0, 0x01, mbpc->comm(), &stat );
1950  MPI_Recv( ids, mySize, MPI_INT, 0, 0x03, mbpc->comm(), &stat );
1951  MPI_Recv( length, mySize, MPI_INT, 0, 0x06, mbpc->comm(), &stat );
1952  sum = 0;
1953 
1954  for( j = 0; j < mySize; j++ )
1955  sum += length[j];
1956 
1957  adjs = (int*)malloc( sizeof( int ) * sum );
1958  if( edge_weights != NULL ) edge_weights = (double*)malloc( sizeof( double ) * sum );
1959  nborProcs = (int*)malloc( sizeof( int ) * sum );
1960  MPI_Recv( adjs, sum, MPI_INT, 0, 0x07, mbpc->comm(), &stat );
1961  MPI_Recv( nborProcs, sum, MPI_INT, 0, 0x08, mbpc->comm(), &stat );
1962  }
1963 
1964  Points = pts;
1965  GlobalIds = ids;
1966  NumPoints = mySize;
1967  NumEdges = length;
1968  NborGlobalId = adjs;
1969  NborProcs = nborProcs;
1970  ObjWeights = obj_weights;
1971  EdgeWeights = edge_weights;
1972  Parts = parts;
1973 
1974  return mySize;
1975 }
1976 
1977 void ZoltanPartitioner::mbFinalizePoints( int npts,
1978  int numExport,
1979  ZOLTAN_ID_PTR exportLocalIDs,
1980  int* exportProcs,
1981  int** assignment )
1982 {
1983  int* MyAssignment;
1984  int i;
1985  int numPts;
1986  MPI_Status stat;
1987  int* recvA;
1988 
1989  /* assign pts to start */
1990 
1991  if( mbpc->proc_config().proc_rank() == 0 )
1992  MyAssignment = (int*)malloc( sizeof( int ) * npts );
1993  else
1994  MyAssignment = (int*)malloc( sizeof( int ) * NumPoints );
1995 
1996  for( i = 0; i < NumPoints; i++ )
1997  MyAssignment[i] = mbpc->proc_config().proc_rank();
1998 
1999  for( i = 0; i < numExport; i++ )
2000  MyAssignment[exportLocalIDs[i]] = exportProcs[i];
2001 
2002  if( mbpc->proc_config().proc_rank() == 0 )
2003  {
2004  /* collect pts */
2005  recvA = MyAssignment + NumPoints;
2006 
2007  for( i = 1; i < (int)mbpc->proc_config().proc_size(); i++ )
2008  {
2009  MPI_Recv( &numPts, 1, MPI_INT, i, 0x04, mbpc->comm(), &stat );
2010  MPI_Recv( recvA, numPts, MPI_INT, i, 0x05, mbpc->comm(), &stat );
2011  recvA += numPts;
2012  }
2013 
2014  *assignment = MyAssignment;
2015  }
2016  else
2017  {
2018  MPI_Send( &NumPoints, 1, MPI_INT, 0, 0x04, mbpc->comm() );
2019  MPI_Send( MyAssignment, NumPoints, MPI_INT, 0, 0x05, mbpc->comm() );
2020  free( MyAssignment );
2021  }
2022 }
2023 
2024 int ZoltanPartitioner::mbGlobalSuccess( int rc )
2025 {
2026  int fail = 0;
2027  unsigned int i;
2028  int* vals = (int*)malloc( mbpc->proc_config().proc_size() * sizeof( int ) );
2029 
2030  MPI_Allgather( &rc, 1, MPI_INT, vals, 1, MPI_INT, mbpc->comm() );
2031 
2032  for( i = 0; i < mbpc->proc_config().proc_size(); i++ )
2033  {
2034  if( vals[i] != ZOLTAN_OK )
2035  {
2036  if( 0 == mbpc->proc_config().proc_rank() )
2037  {
2038  mbShowError( vals[i], "Result on process " );
2039  }
2040  fail = 1;
2041  }
2042  }
2043 
2044  free( vals );
2045  return fail;
2046 }
2047 
2048 void ZoltanPartitioner::mbPrintGlobalResult( const char* s, int begin, int import, int exp, int change )
2049 {
2050  unsigned int i;
2051  int* v1 = (int*)malloc( 4 * sizeof( int ) );
2052  int* v2 = NULL;
2053  int* v;
2054 
2055  v1[0] = begin;
2056  v1[1] = import;
2057  v1[2] = exp;
2058  v1[3] = change;
2059 
2060  if( mbpc->proc_config().proc_rank() == 0 )
2061  {
2062  v2 = (int*)malloc( 4 * mbpc->proc_config().proc_size() * sizeof( int ) );
2063  }
2064 
2065  MPI_Gather( v1, 4, MPI_INT, v2, 4, MPI_INT, 0, mbpc->comm() );
2066 
2067  if( mbpc->proc_config().proc_rank() == 0 )
2068  {
2069  fprintf( stdout, "======%s======\n", s );
2070  for( i = 0, v = v2; i < mbpc->proc_config().proc_size(); i++, v += 4 )
2071  {
2072  fprintf( stdout, "%u: originally had %d, import %d, exp %d, %s\n", i, v[0], v[1], v[2],
2073  v[3] ? "a change of partitioning" : "no change" );
2074  }
2075  fprintf( stdout, "==========================================\n" );
2076  fflush( stdout );
2077 
2078  free( v2 );
2079  }
2080 
2081  free( v1 );
2082 }
2083 
2084 void ZoltanPartitioner::mbShowError( int val, const char* s )
2085 {
2086  if( s ) printf( "%s ", s );
2087 
2088  switch( val )
2089  {
2090  case ZOLTAN_OK:
2091  printf( "%d: SUCCESSFUL\n", mbpc->proc_config().proc_rank() );
2092  break;
2093  case ZOLTAN_WARN:
2094  printf( "%d: WARNING\n", mbpc->proc_config().proc_rank() );
2095  break;
2096  case ZOLTAN_FATAL:
2097  printf( "%d: FATAL ERROR\n", mbpc->proc_config().proc_rank() );
2098  break;
2099  case ZOLTAN_MEMERR:
2100  printf( "%d: MEMORY ALLOCATION ERROR\n", mbpc->proc_config().proc_rank() );
2101  break;
2102  default:
2103  printf( "%d: INVALID RETURN CODE\n", mbpc->proc_config().proc_rank() );
2104  break;
2105  }
2106 }
2107 
2108 /**********************
2109 ** call backs
2110 **********************/
2111 
2112 int mbGetNumberOfAssignedObjects( void* /* userDefinedData */, int* err )
2113 {
2114  *err = 0;
2115  return NumPoints;
2116 }
2117 
2118 void mbGetObjectList( void* /* userDefinedData */,
2119  int /* numGlobalIds */,
2120  int /* numLids */,
2121  ZOLTAN_ID_PTR gids,
2122  ZOLTAN_ID_PTR lids,
2123  int wgt_dim,
2124  float* obj_wgts,
2125  int* err )
2126 {
2127  for( int i = 0; i < NumPoints; i++ )
2128  {
2129  gids[i] = GlobalIds[i];
2130  lids[i] = i;
2131  if( wgt_dim > 0 ) obj_wgts[i] = ObjWeights[i];
2132  }
2133 
2134  *err = 0;
2135 }
2136 
2137 int mbGetObjectSize( void* /* userDefinedData */, int* err )
2138 {
2139  *err = 0;
2140  return 3;
2141 }
2142 
2143 void mbGetObject( void* /* userDefinedData */,
2144  int /* numGlobalIds */,
2145  int /* numLids */,
2146  int numObjs,
2147  ZOLTAN_ID_PTR /* gids */,
2148  ZOLTAN_ID_PTR lids,
2149  int numDim,
2150  double* pts,
2151  int* err )
2152 {
2153  int i, id, id3;
2154  int next = 0;
2155 
2156  if( numDim != 3 )
2157  {
2158  *err = 1;
2159  return;
2160  }
2161 
2162  for( i = 0; i < numObjs; i++ )
2163  {
2164  id = lids[i];
2165 
2166  if( ( id < 0 ) || ( id >= NumPoints ) )
2167  {
2168  *err = 1;
2169  return;
2170  }
2171 
2172  id3 = lids[i] * 3;
2173 
2174  pts[next++] = Points[id3];
2175  pts[next++] = Points[id3 + 1];
2176  pts[next++] = Points[id3 + 2];
2177  }
2178 }
2179 
2180 void mbGetNumberOfEdges( void* /* userDefinedData */,
2181  int /* numGlobalIds */,
2182  int /* numLids */,
2183  int numObjs,
2184  ZOLTAN_ID_PTR /* gids */,
2185  ZOLTAN_ID_PTR lids,
2186  int* numEdges,
2187  int* err )
2188 {
2189  int i, id;
2190  int next = 0;
2191 
2192  for( i = 0; i < numObjs; i++ )
2193  {
2194  id = lids[i];
2195 
2196  if( ( id < 0 ) || ( id >= NumPoints ) )
2197  {
2198  *err = 1;
2199  return;
2200  }
2201 
2202  numEdges[next++] = NumEdges[id];
2203  }
2204 }
2205 
2206 void mbGetEdgeList( void* /* userDefinedData */,
2207  int /* numGlobalIds */,
2208  int /* numLids */,
2209  int numObjs,
2210  ZOLTAN_ID_PTR /* gids */,
2211  ZOLTAN_ID_PTR lids,
2212  int* /* numEdges */,
2213  ZOLTAN_ID_PTR nborGlobalIds,
2214  int* nborProcs,
2215  int wgt_dim,
2216  float* edge_wgts,
2217  int* err )
2218 {
2219  int i, id, idSum, j;
2220  int next = 0;
2221 
2222  for( i = 0; i < numObjs; i++ )
2223  {
2224  id = lids[i];
2225 
2226  if( ( id < 0 ) || ( id >= NumPoints ) )
2227  {
2228  *err = 1;
2229  return;
2230  }
2231 
2232  idSum = 0;
2233 
2234  for( j = 0; j < id; j++ )
2235  idSum += NumEdges[j];
2236 
2237  for( j = 0; j < NumEdges[id]; j++ )
2238  {
2239  nborGlobalIds[next] = NborGlobalId[idSum];
2240  nborProcs[next] = NborProcs[idSum];
2241  if( wgt_dim > 0 ) edge_wgts[next] = EdgeWeights[idSum];
2242  next++;
2243  idSum++;
2244  }
2245  }
2246 }
2247 
2248 void mbGetPart( void* /* userDefinedData */,
2249  int /* numGlobalIds */,
2250  int /* numLids */,
2251  int numObjs,
2252  ZOLTAN_ID_PTR /* gids */,
2253  ZOLTAN_ID_PTR lids,
2254  int* part,
2255  int* err )
2256 {
2257  int i, id;
2258  int next = 0;
2259 
2260  for( i = 0; i < numObjs; i++ )
2261  {
2262  id = lids[i];
2263 
2264  if( ( id < 0 ) || ( id >= NumPoints ) )
2265  {
2266  *err = 1;
2267  return;
2268  }
2269 
2270  part[next++] = Parts[id];
2271  }
2272 }
2273 
2274 // new methods for partition in parallel, used by migrate in iMOAB
2275 ErrorCode ZoltanPartitioner::partition_owned_cells( Range& primary,
2276  std::multimap< int, int >& extraGraphEdges,
2277  std::map< int, int > procs,
2278  int& numNewPartitions,
2279  std::map< int, Range >& distribution,
2280  int met )
2281 {
2282  // start copy
2283  MeshTopoUtil mtu( mbImpl );
2284  Range adjs;
2285 
2286  std::vector< int > adjacencies;
2287  std::vector< int > ids;
2288  ids.resize( primary.size() );
2289  std::vector< int > length;
2290  std::vector< double > coords;
2291  std::vector< int > nbor_proc;
2292  // can use a fixed-size array 'cuz the number of lower-dimensional neighbors is limited
2293  // by MBCN
2294  int neighbors[5 * MAX_SUB_ENTITIES];
2295  int neib_proc[5 * MAX_SUB_ENTITIES];
2296  double avg_position[3];
2297  int moab_id;
2298  int primaryDim = mbImpl->dimension_from_handle( *primary.rbegin() );
2299  // get the global id tag handle
2300  Tag gid = mbImpl->globalId_tag();
2301 
2302  ErrorCode rval = mbImpl->tag_get_data( gid, primary, &ids[0] );MB_CHK_ERR( rval );
2303 
2304  // mbpc is member in base class, PartitionerBase
2305  int rank = mbpc->rank(); // current rank , will be put on regular neighbors
2306  int i = 0;
2307  for( Range::iterator rit = primary.begin(); rit != primary.end(); ++rit, i++ )
2308  {
2309  EntityHandle cell = *rit;
2310  // get bridge adjacencies for each cell
2311  if( 1 == met )
2312  {
2313  adjs.clear();
2314  rval = mtu.get_bridge_adjacencies( cell, ( primaryDim > 0 ? primaryDim - 1 : 3 ), primaryDim, adjs );MB_CHK_ERR( rval );
2315 
2316  // get the graph vertex ids of those
2317  if( !adjs.empty() )
2318  {
2319  assert( adjs.size() < 5 * MAX_SUB_ENTITIES );
2320  rval = mbImpl->tag_get_data( gid, adjs, neighbors );MB_CHK_ERR( rval );
2321  }
2322  // if adjacent to neighbor partitions, add to the list
2323  int size_adjs = (int)adjs.size();
2324  moab_id = ids[i];
2325  for( int k = 0; k < size_adjs; k++ )
2326  neib_proc[k] = rank; // current rank
2327  if( extraGraphEdges.find( moab_id ) != extraGraphEdges.end() )
2328  {
2329  // it means that the current cell is adjacent to a cell in another partition ; maybe
2330  // a few
2331  std::pair< std::multimap< int, int >::iterator, std::multimap< int, int >::iterator > ret;
2332  ret = extraGraphEdges.equal_range( moab_id );
2333  for( std::multimap< int, int >::iterator it = ret.first; it != ret.second; ++it )
2334  {
2335  int otherID = it->second;
2336  neighbors[size_adjs] = otherID; // the id of the other cell, across partition
2337  neib_proc[size_adjs] = procs[otherID]; // this is how we built this map, the cell id maps to
2338  // what proc it came from
2339  size_adjs++;
2340  }
2341  }
2342  // copy those into adjacencies vector
2343  length.push_back( size_adjs );
2344  std::copy( neighbors, neighbors + size_adjs, std::back_inserter( adjacencies ) );
2345  std::copy( neib_proc, neib_proc + size_adjs, std::back_inserter( nbor_proc ) );
2346  }
2347  else if( 2 == met )
2348  {
2349  if( TYPE_FROM_HANDLE( cell ) == MBVERTEX )
2350  {
2351  rval = mbImpl->get_coords( &cell, 1, avg_position );MB_CHK_ERR( rval );
2352  }
2353  else
2354  {
2355  rval = mtu.get_average_position( cell, avg_position );MB_CHK_ERR( rval );
2356  }
2357  std::copy( avg_position, avg_position + 3, std::back_inserter( coords ) );
2358  }
2359  }
2360 
2361 #ifdef VERBOSE
2362  std::stringstream ff2;
2363  ff2 << "zoltanInput_" << mbpc->rank() << ".txt";
2364  std::ofstream ofs;
2365  ofs.open( ff2.str().c_str(), std::ofstream::out );
2366  ofs << "Length vector: " << std::endl;
2367  std::copy( length.begin(), length.end(), std::ostream_iterator< int >( ofs, ", " ) );
2368  ofs << std::endl;
2369  ofs << "Adjacencies vector: " << std::endl;
2370  std::copy( adjacencies.begin(), adjacencies.end(), std::ostream_iterator< int >( ofs, ", " ) );
2371  ofs << std::endl;
2372  ofs << "Moab_ids vector: " << std::endl;
2373  std::copy( ids.begin(), ids.end(), std::ostream_iterator< int >( ofs, ", " ) );
2374  ofs << std::endl;
2375  ofs << "Coords vector: " << std::endl;
2376  std::copy( coords.begin(), coords.end(), std::ostream_iterator< double >( ofs, ", " ) );
2377  ofs << std::endl;
2378  ofs.close();
2379 #endif
2380 
2381  // these are static var in this file, and used in the callbacks
2382  Points = NULL;
2383  if( 1 != met ) Points = &coords[0];
2384  GlobalIds = &ids[0];
2385  NumPoints = (int)ids.size();
2386  NumEdges = &length[0];
2387  NborGlobalId = &adjacencies[0];
2388  NborProcs = &nbor_proc[0];
2389  ObjWeights = NULL;
2390  EdgeWeights = NULL;
2391  Parts = NULL;
2392 
2393  float version;
2394  if( mbpc->rank() == 0 ) std::cout << "Initializing zoltan..." << std::endl;
2395 
2396  Zoltan_Initialize( argcArg, argvArg, &version );
2397 
2398  // Create Zoltan object. This calls Zoltan_Create.
2399  if( NULL == myZZ ) myZZ = new Zoltan( mbpc->comm() );
2400 
2401  // set # requested partitions
2402  char buff[10];
2403  sprintf( buff, "%d", numNewPartitions );
2404  int retval = myZZ->Set_Param( "NUM_GLOBAL_PARTITIONS", buff );
2405  if( ZOLTAN_OK != retval ) return MB_FAILURE;
2406 
2407  // request parts assignment
2408  retval = myZZ->Set_Param( "RETURN_LISTS", "PARTS" );
2409  if( ZOLTAN_OK != retval ) return MB_FAILURE;
2410 
2411  myZZ->Set_Num_Obj_Fn( mbGetNumberOfAssignedObjects, NULL );
2412  myZZ->Set_Obj_List_Fn( mbGetObjectList, NULL );
2413  // due to a bug in zoltan, if method is graph partitioning, do not pass coordinates!!
2414  if( 2 == met )
2415  {
2416  myZZ->Set_Num_Geom_Fn( mbGetObjectSize, NULL );
2417  myZZ->Set_Geom_Multi_Fn( mbGetObject, NULL );
2418  SetRCB_Parameters(); // geometry
2419  }
2420  else if( 1 == met )
2421  {
2422  myZZ->Set_Num_Edges_Multi_Fn( mbGetNumberOfEdges, NULL );
2423  myZZ->Set_Edge_List_Multi_Fn( mbGetEdgeList, NULL );
2424  SetHypergraph_Parameters( "auto" );
2425  }
2426 
2427  // Perform the load balancing partitioning
2428 
2429  int changes;
2430  int numGidEntries;
2431  int numLidEntries;
2432  int num_import;
2433  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
2434  int* import_procs;
2435  int* import_to_part;
2436  int num_export;
2437  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
2438  int *assign_procs, *assign_parts;
2439 
2440  if( mbpc->rank() == 0 )
2441  std::cout << "Computing partition using method (1-graph, 2-geom):" << met << " for " << numNewPartitions
2442  << " parts..." << std::endl;
2443 
2444 #ifndef NDEBUG
2445 #if 0
2446  static int counter=0; // it may be possible to call function multiple times in a simulation
2447  // give a way to not overwrite the files
2448  // it should work only with a modified version of Zoltan
2449  std::stringstream basename;
2450  if (1==met)
2451  {
2452  basename << "phg_" << counter++;
2453  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 0, 1, 0);
2454  }
2455  else if (2==met)
2456  {
2457  basename << "rcb_" << counter++;
2458  Zoltan_Generate_Files(myZZ->Get_C_Handle(), (char*)(basename.str().c_str()), 1, 1, 0, 0);
2459  }
2460 #endif
2461 #endif
2462  retval = myZZ->LB_Partition( changes, numGidEntries, numLidEntries, num_import, import_global_ids, import_local_ids,
2463  import_procs, import_to_part, num_export, export_global_ids, export_local_ids,
2464  assign_procs, assign_parts );
2465  if( ZOLTAN_OK != retval ) return MB_FAILURE;
2466 
2467 #ifdef VERBOSE
2468  std::stringstream ff3;
2469  ff3 << "zoltanOutput_" << mbpc->rank() << ".txt";
2470  std::ofstream ofs3;
2471  ofs3.open( ff3.str().c_str(), std::ofstream::out );
2472  ofs3 << " export elements on rank " << rank << " \n";
2473  ofs3 << "\t index \t gb_id \t local \t proc \t part \n";
2474  for( int k = 0; k < num_export; k++ )
2475  {
2476  ofs3 << "\t" << k << "\t" << export_global_ids[k] << "\t" << export_local_ids[k] << "\t" << assign_procs[k]
2477  << "\t" << assign_parts[k] << "\n";
2478  }
2479  ofs3.close();
2480 #endif
2481 
2482  // basically each local cell is assigned to a part
2483 
2484  assert( num_export == (int)primary.size() );
2485  for( i = 0; i < num_export; i++ )
2486  {
2487  EntityHandle cell = primary[export_local_ids[i]];
2488  distribution[assign_parts[i]].insert( cell );
2489  }
2490 
2491  Zoltan::LB_Free_Part( &import_global_ids, &import_local_ids, &import_procs, &import_to_part );
2492  Zoltan::LB_Free_Part( &export_global_ids, &export_local_ids, &assign_procs, &assign_parts );
2493 
2494  delete myZZ;
2495  myZZ = NULL;
2496 
2497  // clear arrays that were resized locally, to free up local memory
2498 
2499  std::vector< int >().swap( adjacencies );
2500  std::vector< int >().swap( ids );
2501  std::vector< int >().swap( length );
2502  std::vector< int >().swap( nbor_proc );
2503  std::vector< double >().swap( coords );
2504 
2505  return MB_SUCCESS;
2506 }