Actual source code: Partitioner.hh

  1: #ifndef included_ALE_Partitioner_hh
  2: #define included_ALE_Partitioner_hh

  4: #ifndef  included_ALE_Completion_hh
  5: #include <Completion.hh>
  6: #endif

  8: #ifdef PETSC_HAVE_ZOLTAN
  9: #include <zoltan.h>

 12:   // Inputs

 18:   int getNumVertices_Zoltan(void *, int *);

 20:   void getLocalElements_Zoltan(void *, int, int, ZOLTAN_ID_PTR, ZOLTAN_ID_PTR, int, float *, int *);

 22:   void getHgSizes_Zoltan(void *, int *, int *, int *, int *);

 24:   void getHg_Zoltan(void *, int, int, int, int, ZOLTAN_ID_PTR, int *, ZOLTAN_ID_PTR, int *);
 25: }

 27: #endif

 29: #ifdef PETSC_HAVE_CHACO
 30: #ifdef PETSC_HAVE_UNISTD_H
 31: #include <unistd.h>
 32: #endif
 33: /* Chaco does not have an include file */
 36:                        float *ewgts, float *x, float *y, float *z, char *outassignname,
 37:                        char *outfilename, short *assignment, int architecture, int ndims_tot,
 38:                        int mesh_dims[3], double *goal, int global_method, int local_method,
 39:                        int rqi_flag, int vmax, int ndims, double eigtol, long seed);

 42: }
 43: #endif
 44: #ifdef PETSC_HAVE_PARMETIS
 46:   #include <parmetis.h>
 48: }
 49: #endif
 50: #ifdef PETSC_HAVE_HMETIS
 53: }
 54: #endif

 56: namespace ALE {
 57: #if 1
 58: #ifdef PETSC_HAVE_CHACO
 59:   namespace Chaco {
 60:     template<typename Alloc_ = malloc_allocator<short int> >
 61:     class Partitioner {
 62:     public:
 63:       typedef short int part_type;
 64:       typedef Alloc_    alloc_type;
 65:       enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
 66:     protected:
 67:       const int  logSize;
 68:       char      *msgLog;
 69:       int        fd_stdout, fd_pipe[2];
 70:       alloc_type _allocator;
 71:     public:
 72:       Partitioner(): logSize(10000) {};
 73:       ~Partitioner() {};
 74:     protected:
 75:       // Chaco outputs to stdout. We redirect this to a buffer.
 76:       // TODO: check error codes for UNIX calls
 77:       void startStdoutRedirect() {
 78: #ifdef PETSC_HAVE_UNISTD_H
 79:         this->fd_stdout = dup(1);
 80:         pipe(this->fd_pipe);
 81:         close(1);
 82:         dup2(this->fd_pipe[1], 1);
 83: #endif
 84:       };
 85:       void stopStdoutRedirect() {
 86: #ifdef PETSC_HAVE_UNISTD_H
 87:         int count;

 89:         fflush(stdout);
 90:         this->msgLog = new char[this->logSize];
 91:         count = read(this->fd_pipe[0], this->msgLog, (this->logSize-1)*sizeof(char));
 92:         if (count < 0) count = 0;
 93:         this->msgLog[count] = 0;
 94:         close(1);
 95:         dup2(this->fd_stdout, 1);
 96:         close(this->fd_stdout);
 97:         close(this->fd_pipe[0]);
 98:         close(this->fd_pipe[1]);
 99:         //std::cout << this->msgLog << std::endl;
100:         delete [] this->msgLog;
101: #endif
102:       };
103:     public:
104:       static bool zeroBase() {return false;}
105:       // This method returns the partition section mapping sieve points (here cells) to partitions
106:       //   start:     start of edge list for each vertex
107:       //   adjacency: adj[start[v]] is edge list data for vertex v
108:       //   partition: this section is over the partitions and takes points as values
109:       // TODO: Read global and local methods from options
110:       template<typename Section, typename MeshManager>
111:       void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
112:         FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
113:         int nvtxs = numVertices;                /* number of vertices in full graph */
114:         int *vwgts = NULL;                      /* weights for all vertices */
115:         float *ewgts = NULL;                    /* weights for all edges */
116:         float *x = NULL, *y = NULL, *z = NULL;  /* coordinates for inertial method */
117:         char *outassignname = NULL;             /*  name of assignment output file */
118:         char *outfilename = NULL;               /* output file name */
119:         int architecture = 1;                   /* 0 => hypercube, d => d-dimensional mesh */
120:         int ndims_tot = 0;                      /* total number of cube dimensions to divide */
121:         int mesh_dims[3];                       /* dimensions of mesh of processors */
122:         double *goal = NULL;                    /* desired set sizes for each set */
123:         int global_method = 1;                  /* global partitioning algorithm */
124:         int local_method = 1;                   /* local partitioning algorithm */
125:         int rqi_flag = 0;                       /* should I use RQI/Symmlq eigensolver? */
126:         int vmax = 200;                         /* how many vertices to coarsen down to? */
127:         int ndims = 1;                          /* number of eigenvectors (2^d sets) */
128:         double eigtol = 0.001;                  /* tolerance on eigenvectors */
129:         long seed = 123636512;                  /* for random graph mutations */
130:         int maxSize = 0;

132:             if (global_method == INERTIAL_METHOD) {manager.createCellCoordinates(nvtxs, &x, &y, &z);}
133:         mesh_dims[0] = partition->commSize(); mesh_dims[1] = 1; mesh_dims[2] = 1;
134:         part_type *assignment = this->_allocator.allocate(nvtxs);
135:         for(int i = 0; i < nvtxs; ++i) {this->_allocator.construct(assignment+i, 0);}

137:         this->startStdoutRedirect();
138:         interface(nvtxs, (int *) start, (int *) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
139:                   assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
140:                   vmax, ndims, eigtol, seed);
141:         this->stopStdoutRedirect();

143:         for(int v = 0; v < nvtxs; ++v) {partition->addFiberDimension(assignment[v], 1);}
144:         partition->allocatePoint();
145:         for(int p = 0; p < partition->commSize(); ++p) {
146:           maxSize = std::max(maxSize, partition->getFiberDimension(p));
147:         }
148:         typename Section::value_type *values = new typename Section::value_type[maxSize];

150:         for(int p = 0; p < partition->commSize(); ++p) {
151:           int k = 0;

153:           for(int v = 0; v < nvtxs; ++v) {
154:             if (assignment[v] == p) values[k++] = manager.getCell(v);
155:           }
156:           if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
157:           partition->updatePoint(p, values);
158:         }
159:         delete [] values;

161:             if (global_method == INERTIAL_METHOD) {manager.destroyCellCoordinates(nvtxs, &x, &y, &z);}
162:         for(int i = 0; i < nvtxs; ++i) {this->_allocator.destroy(assignment+i);}
163:         this->_allocator.deallocate(assignment, nvtxs);
164:       };
165:     };
166:   }
167: #endif
168: #ifdef PETSC_HAVE_PARMETIS
169:   namespace ParMetis {
170:     template<typename Alloc_ = malloc_allocator<int> >
171:     class Partitioner {
172:     public:
173:       typedef int    part_type;
174:       typedef Alloc_ alloc_type;
175:     protected:
176:       alloc_type _allocator;
177:     public:
178:       static bool zeroBase() {return true;}
179:       // This method returns the partition section mapping sieve points (here cells) to partitions
180:       //   start:     start of edge list for each vertex
181:       //   adjacency: adj[start[v]] is edge list data for vertex v
182:       //   partition: this section is over the partitions and takes points as values
183:       // TODO: Read parameters from options
184:       template<typename Section, typename MeshManager>
185:       void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
186:         //static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
187:         int    nvtxs      = numVertices; // The number of vertices in full graph
188:         int   *vtxdist;                  // Distribution of vertices across processes
189:         int   *xadj       = const_cast<int*>(start);       // Start of edge list for each vertex
190:         int   *adjncy     = const_cast<int*>(adjacency);   // Edge lists for all vertices
191:         int   *vwgt       = NULL;        // Vertex weights
192:         int   *adjwgt     = NULL;        // Edge weights
193:         int    wgtflag    = 0;           // Indicates which weights are present
194:         int    numflag    = 0;           // Indicates initial offset (0 or 1)
195:         int    ncon       = 1;           // The number of weights per vertex
196:         int    nparts     = partition->commSize(); // The number of partitions
197:         float *tpwgts;                   // The fraction of vertex weights assigned to each partition
198:         float *ubvec;                    // The balance intolerance for vertex weights
199:         int    options[5];               // Options
200:         int    maxSize    = 0;
201:         // Outputs
202:         int        edgeCut;              // The number of edges cut by the partition
203:         part_type *assignment;

205:         options[0] = 0; // Use all defaults
206:         // Calculate vertex distribution
207:         //   Not sure this still works in parallel
208:         vtxdist    = new int[nparts+1];
209:         vtxdist[0] = 0;
210:         MPI_Allgather(&nvtxs, 1, MPI_INT, &vtxdist[1], 1, MPI_INT, partition->comm());
211:         for(int p = 2; p <= nparts; ++p) {
212:           vtxdist[p] += vtxdist[p-1];
213:         }
214:         // Calculate weights
215:         tpwgts     = new float[ncon*nparts];
216:         for(int p = 0; p < nparts; ++p) {
217:           tpwgts[p] = 1.0/nparts;
218:         }
219:         ubvec      = new float[ncon];
220:         ubvec[0]   = 1.05;

222:         assignment = this->_allocator.allocate(nvtxs);
223:         for(int i = 0; i < nvtxs; ++i) {this->_allocator.construct(assignment+i, 0);}

225:         if (partition->commSize() == 1) {
226:           PetscMemzero(assignment, nvtxs * sizeof(part_type));
227:         } else {
228:           if (partition->debug() && nvtxs) {
229:             for(int p = 0; p <= nvtxs; ++p) {
230:               std::cout << "["<<partition->commRank()<<"]xadj["<<p<<"] = " << xadj[p] << std::endl;
231:             }
232:             for(int i = 0; i < xadj[nvtxs]; ++i) {
233:               std::cout << "["<<partition->commRank()<<"]adjncy["<<i<<"] = " << adjncy[i] << std::endl;
234:             }
235:           }
236:           if (vtxdist[1] == vtxdist[nparts]) {
237:             if (partition->commRank() == 0) {
238:               METIS_PartGraphKway(&nvtxs, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &nparts, options, &edgeCut, assignment);
239:               if (partition->debug()) {std::cout << "Metis: edgecut is " << edgeCut << std::endl;}
240:             }
241:           } else {
242:             MPI_Comm comm = partition->comm();

244:             ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
245:             if (partition->debug()) {std::cout << "ParMetis: edgecut is " << edgeCut << std::endl;}
246:           }
247:         }
248:         delete [] vtxdist;
249:         delete [] tpwgts;
250:         delete [] ubvec;

252:         for(int v = 0; v < nvtxs; ++v) {partition->addFiberDimension(assignment[v], 1);}
253:         partition->allocatePoint();
254:         for(int p = 0; p < partition->commSize(); ++p) {
255:           maxSize = std::max(maxSize, partition->getFiberDimension(p));
256:         }
257:         typename Section::value_type *values = new typename Section::value_type[maxSize];

259:         for(int p = 0; p < partition->commSize(); ++p) {
260:           int k = 0;

262:           for(int v = 0; v < nvtxs; ++v) {
263:             if (assignment[v] == p) values[k++] = manager.getCell(v);
264:           }
265:           if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
266:           partition->updatePoint(p, values);
267:         }
268:         delete [] values;

270:         for(int i = 0; i < nvtxs; ++i) {this->_allocator.destroy(assignment+i);}
271:         this->_allocator.deallocate(assignment, nvtxs);
272:       };
273:     };
274:   };
275: #endif
276:   namespace Simple {
277:     template<typename Alloc_ = malloc_allocator<short int> >
278:     class Partitioner {
279:     public:
280:       typedef int    part_type;
281:       typedef Alloc_ alloc_type;
282:     protected:
283:       alloc_type _allocator;
284:     public:
285:       Partitioner() {};
286:       ~Partitioner() {};
287:     public:
288:       static bool zeroBase() {return true;}
289:       template<typename Section, typename MeshManager>
290:       void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
291:         const int numProcs = partition->commSize();
292:         const int rank     = partition->commRank();
293:         int       maxSize  = 0;

295:         for(int p = 0; p < numProcs; ++p) {
296:           partition->setFiberDimension(p, numVertices/numProcs + ((numVertices % numProcs) > rank));
297:           maxSize = std::max(maxSize, partition->getFiberDimension(p));
298:         }
299:         partition->allocatePoint();
300:         typename Section::value_type *values = new typename Section::value_type[maxSize];

302:         for(int p = 0; p < partition->commSize(); ++p) {
303:           const int start = p*(numVertices/numProcs)     + p*((numVertices % numProcs) > p+1);
304:           const int end   = (p+1)*(numVertices/numProcs) + (p+1)*((numVertices % numProcs) > p+2);
305:           int       k     = 0;

307:           for(int v = start; v < end; ++v, ++k) {
308:             values[k] = manager.getCell(v);
309:           }
310:           if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
311:           partition->updatePoint(p, values);
312:         }
313:         delete [] values;
314:       };
315:     };
316:   }
317: #ifdef PETSC_HAVE_CHACO
318:   template<typename GraphPartitioner = ALE::Chaco::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
319: #elif defined(PETSC_HAVE_PARMETIS)
320:   template<typename GraphPartitioner = ALE::ParMetis::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
321: #else
322:   template<typename GraphPartitioner = ALE::Simple::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
323: #endif
324:   class Partitioner {
325:   public:
326:     typedef Alloc_                               alloc_type;
327:     typedef GraphPartitioner                     graph_partitioner_type;
328:     typedef typename GraphPartitioner::part_type part_type;
329:     template<typename Mesh>
330:     class MeshManager {
331:     public:
332:       typedef typename Mesh::point_type point_type;
333:     protected:
334:       const Obj<Mesh>& mesh;
335:       bool             simpleCellNumbering;
336:       point_type      *cells;
337:     protected:
338:       void createCells(const int height) {
339:         const Obj<typename Mesh::label_sequence>&     mcells   = mesh->heightStratum(height);
340:         const typename Mesh::label_sequence::iterator cEnd     = mcells->end();
341:         const int                                     numCells = mcells->size();
342:         int                                           c        = 0;

344:         this->cells               = NULL;
345:         this->simpleCellNumbering = true;
346:         for(typename Mesh::label_sequence::iterator c_iter = mcells->begin(); c_iter != cEnd; ++c_iter, ++c) {
347:           if (*c_iter != c) {
348:             this->simpleCellNumbering = false;
349:             break;
350:           }
351:         }
352:         if (!this->simpleCellNumbering) {
353:           this->cells = new point_type[numCells];
354:           c           = 0;
355:           for(typename Mesh::label_sequence::iterator c_iter = mcells->begin(); c_iter != cEnd; ++c_iter, ++c) {
356:             this->cells[c] = *c_iter;
357:           }
358:         }
359:       };
360:     public:
361:       MeshManager(const Obj<Mesh>& mesh, const int height = 0): mesh(mesh) {
362:         this->createCells(height);
363:       };
364:       ~MeshManager() {
365:         if (this->cells) {delete [] this->cells;}
366:       };
367:     public:
368:       template<typename Float>
369:       void createCellCoordinates(const int numVertices, Float *X[], Float *Y[], Float *Z[]) const {
370:         const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
371:         const int dim = mesh->getDimension();
372:         typedef typename alloc_type::template rebind<Float>::other float_alloc_type;
373:         Float *x = float_alloc_type().allocate(numVertices*3);
374:         for(int i = 0; i < numVertices*3; ++i) {float_alloc_type().construct(x+i, 0.0);}
375:         Float *y = x+numVertices;
376:         Float *z = y+numVertices;
377:         Float *vCoords[3];

379:         vCoords[0] = x; vCoords[1] = y; vCoords[2] = z;
380:         const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
381:         const int corners = mesh->size(coordinates, *(cells->begin()))/dim;
382:         int       c       = 0;

384:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter !=cells->end(); ++c_iter, ++c) {
385:           const double *coords = mesh->restrictClosure(coordinates, *c_iter);

387:           for(int d = 0; d < dim; ++d) {
388:             vCoords[d][c] = 0.0;
389:           }
390:           for(int v = 0; v < corners; ++v) {
391:             for(int d = 0; d < dim; ++d) {
392:               vCoords[d][c] += coords[v*dim+d];
393:             }
394:           }
395:           for(int d = 0; d < dim; ++d) {
396:             vCoords[d][c] /= corners;
397:           }
398:         }
399:         *X = x;
400:         *Y = y;
401:         *Z = z;
402:       };
403:       template<typename Float>
404:       void destroyCellCoordinates(const int numVertices, Float *X[], Float *Y[], Float *Z[]) const {
405:         typedef typename alloc_type::template rebind<Float>::other float_alloc_type;
406:         Float *x = *X;

408:         for(int i = 0; i < numVertices*3; ++i) {float_alloc_type().destroy(x+i);}
409:         float_alloc_type().deallocate(x, numVertices*3);
410:       };
411:       point_type getCell(const int cellNumber) const {
412:         if (this->simpleCellNumbering) {
413:           return cellNumber;
414:         }
415:         return this->cells[cellNumber];
416:       };
417:     };
418:     template<typename Sieve>
419:     class OffsetVisitor {
420:       const Sieve& sieve;
421:       const Sieve& overlapSieve;
422:       int         *offsets;
423:     public:
424:       OffsetVisitor(const Sieve& s, const Sieve& ovS, int off[]) : sieve(s), overlapSieve(ovS), offsets(off) {};
425:       void visitPoint(const typename Sieve::point_type& point) {};
426:       void visitArrow(const typename Sieve::arrow_type& arrow) {
427:         const typename Sieve::point_type cell   = arrow.target;
428:         const typename Sieve::point_type face   = arrow.source;
429:         const int                        size   = this->sieve.getSupportSize(face);
430:         const int                        ovSize = this->overlapSieve.getSupportSize(face);

432:         if (size == 2) {
433:           offsets[cell+1]++;
434:         } else if ((size == 1) && (ovSize == 1)) {
435:           offsets[cell+1]++;
436:         }
437:       };
438:     };
439:     template<typename Sieve>
440:     class AdjVisitor {
441:     protected:
442:       typename Sieve::point_type cell;
443:       int                       *adjacency;
444:       const int                  cellOffset;
445:       int                        offset;
446:     public:
447:       AdjVisitor(int adj[], const bool zeroBase) : adjacency(adj), cellOffset(zeroBase ? 0 : 1), offset(0) {};
448:       void visitPoint(const typename Sieve::point_type& point) {};
449:       void visitArrow(const typename Sieve::arrow_type& arrow) {
450:         const int neighbor = arrow.target;

452:         if (neighbor != this->cell) {
453:           //std::cout << "Adding dual edge from " << cell << " to " << neighbor << std::endl;
454:           this->adjacency[this->offset++] = neighbor + this->cellOffset;
455:         }
456:       };
457:     public:
458:       void setCell(const typename Sieve::point_type cell) {this->cell = cell;};
459:       int  getOffset() {return this->offset;}
460:     };
461:     template<typename Sieve>
462:     class MeetVisitor {
463:     public:
464:       typedef std::set<typename Sieve::point_type> neighbors_type;
465:     protected:
466:       const Sieve& sieve;
467:       const int numCells;
468:       const int faceVertices;
469:       typename Sieve::point_type newCell;
470:       neighbors_type *neighborCells;
471:       typename ISieveVisitor::PointRetriever<Sieve> *pR;
472:       typename Sieve::point_type cell;
473:       std::map<typename Sieve::point_type, typename Sieve::point_type> newCells;
474:     public:
475:       MeetVisitor(const Sieve& s, const int n, const int fV) : sieve(s), numCells(n), faceVertices(fV), newCell(n) {
476:         this->neighborCells = new std::set<typename Sieve::point_type>[numCells];
477:         this->pR            = new typename ISieveVisitor::PointRetriever<Sieve>(this->sieve.getMaxConeSize());
478:       };
479:       ~MeetVisitor() {delete [] this->neighborCells; delete this->pR;};
480:       void visitArrow(const typename Sieve::arrow_type& arrow) {};
481:       void visitPoint(const typename Sieve::point_type& point) {
482:         const typename Sieve::point_type& neighbor = point;

484:         if (this->cell == neighbor) return;
485:         this->pR->clear();
486:         this->sieve.meet(this->cell, neighbor, *this->pR);
487:         if (this->pR->getSize() == (size_t) this->faceVertices) {
488:           if ((this->cell < numCells) && (neighbor < numCells)) {
489:             this->neighborCells[this->cell].insert(neighbor);
490:           } else {
491:             typename Sieve::point_type e = this->cell, n = neighbor;

493:             if (this->cell >= numCells) {
494:               if (this->newCells.find(cell) == this->newCells.end()) this->newCells[cell] = --newCell;
495:               e = this->newCells[cell];
496:             }
497:             if (neighbor >= numCells) {
498:               if (this->newCells.find(neighbor) == this->newCells.end()) this->newCells[neighbor] = --newCell;
499:               n = this->newCells[neighbor];
500:             }
501:             this->neighborCells[e].insert(n);
502:           }
503:         }
504:       };
505:     public:
506:       void setCell(const typename Sieve::point_type& c) {this->cell = c;};
507:       const neighbors_type *getNeighbors() {return this->neighborCells;};
508:     };
509:   public: // Creating overlaps
510:     // Create a partition point overlap for distribution
511:     //   This is the default overlap which comes from distributing a serial mesh on process 0
512:     template<typename SendOverlap, typename RecvOverlap>
513:     static void createDistributionPartOverlap(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
514:       const int rank = sendOverlap->commRank();
515:       const int size = sendOverlap->commSize();

517:       if (rank == 0) {
518:         for(int p = 1; p < size; p++) {
519:           // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
520:           sendOverlap->addCone(p, p, p);
521:         }
522:       }
523:       if (rank != 0) {
524:         // The arrow is from remote partition point rank (color) on rank 0 (source) to local partition point rank (target)
525:         recvOverlap->addCone(0, rank, rank);
526:       }
527:     };
528:     // Create a mesh point overlap for distribution
529:     //   A local numbering is created for the remote points
530:     //   This is the default overlap which comes from distributing a serial mesh on process 0
531:     template<typename Section, typename RecvPartOverlap, typename Renumbering, typename SendOverlap, typename RecvOverlap>
532:     static void createDistributionMeshOverlap(const Obj<Section>& partition, const Obj<RecvPartOverlap>& recvPartOverlap, Renumbering& renumbering, const Obj<Section>& overlapPartition, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
533:       const typename Section::chart_type& chart = partition->getChart();

535:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
536:         if (*p_iter == sendOverlap->commRank()) continue;
537:         const typename Section::value_type *points     = partition->restrictPoint(*p_iter);
538:         const int                           numPoints  = partition->getFiberDimension(*p_iter);
539: 
540:         for(int i = 0; i < numPoints; ++i) {
541:           // Notice here that we do not know the local renumbering (but we do not use it)
542:           sendOverlap->addArrow(points[i], *p_iter, points[i]);
543:         }
544:       }
545:       if (sendOverlap->debug()) {sendOverlap->view("Send mesh overlap");}
546:       const Obj<typename RecvPartOverlap::traits::baseSequence> rPoints    = recvPartOverlap->base();

548:       for(typename RecvPartOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
549:         const Obj<typename RecvPartOverlap::coneSequence>& ranks           = recvPartOverlap->cone(*p_iter);
550:         //const typename Section::point_type&                localPartPoint  = *p_iter;
551:         const typename Section::point_type                 rank            = *ranks->begin();
552:         const typename Section::point_type&                remotePartPoint = ranks->begin().color();
553:         const typename Section::value_type                *points          = overlapPartition->restrictPoint(remotePartPoint);
554:         const int                                          numPoints       = overlapPartition->getFiberDimension(remotePartPoint);

556:         for(int i = 0; i < numPoints; ++i) {
557:           recvOverlap->addArrow(rank, renumbering[points[i]], points[i]);
558:         }
559:       }
560:       if (recvOverlap->debug()) {recvOverlap->view("Receive mesh overlap");}
561:     };
562:     // Create a partition point overlap from a partition
563:     //   The intention is to create an overlap which enables exchange of redistribution information
564:     template<typename Section, typename SendOverlap, typename RecvOverlap>
565:     static void createPartitionPartOverlap(const Obj<Section>& partition, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
566:       const typename Section::chart_type& chart      = partition->getChart();
567:       const int                           rank       = partition->commRank();
568:       const int                           size       = partition->commSize();
569:       int                                *adj        = new int[size];
570:       int                                *recvCounts = new int[size];
571:       int                                 numNeighbors;

573:       for(int p = 0; p < size; ++p) {
574:         adj[p]        = 0;
575:         recvCounts[p] = 1;
576:       }
577:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart->end(); ++p_iter) {
578:         const typename Section::value_type& p = partition->restrictPoint(*p_iter)[0];
579:         // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
580:         sendOverlap->addCone(p, p, p);
581:         adj[p] = 1;
582:       }
583:       MPI_Reduce_Scatter(adj, &numNeighbors, recvCounts, size, MPI_INT, MPI_SUM, partition->comm());
584:       MPI_Request *recvRequests = new MPI_Request[numNeighbors];
585:       int          dummy        = 0;

587:       // TODO: Get a unique tag
588:       for(int n = 0; n < numNeighbors; ++n) {
589:         MPI_Irecv(&dummy, 1, MPI_INT, MPI_ANY_SOURCE, 1, partition->comm(), &recvRequests[n]);
590:       }
591:       const Obj<typename SendOverlap::traits::baseSequence>      ranks        = sendOverlap->base();
592:       const typename SendOverlap::traits::baseSequence::iterator rEnd         = ranks->end();
593:       MPI_Request                                               *sendRequests = new MPI_Request[ranks->size()];
594:       int                                                        s            = 0;

596:       for(typename SendOverlap::traits::baseSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter, ++s) {
597:         MPI_Isend(&dummy, 1, MPI_INT, *r_iter, 1, partition->comm(), &sendRequests[s]);
598:       }
599:       for(int n = 0; n < numNeighbors; ++n) {
600:         MPI_Status status;
601:         int        idx;

603:         MPI_Waitany(numNeighbors, recvRequests, &idx, &status);
604:         // The arrow is from remote partition point rank (color) on rank p (source) to local partition point rank (target)
605:         recvOverlap->addCone(status.MPI_SOURCE, rank, rank);
606:       }
607:       MPI_Waitall(ranks->size(), sendRequests, MPI_STATUSES_IGNORE);
608:       delete [] sendRequests;
609:       delete [] recvRequests;
610:       delete [] adj;
611:       delete [] recvCounts;
612:     };
613:   public: // Building CSR meshes
614:     // This produces the dual graph (each cell is a vertex and each face is an edge)
615:     //   numbering:   A contiguous numbering of the cells (not yet included)
616:     //   numVertices: The number of vertices in the graph (cells in the mesh)
617:     //   adjacency:   The vertices adjacent to each vertex (cells adjacent to each mesh cell)
618:     // - We allow an exception to contiguous numbering.
619:     //   If the cell id > numElements, we assign a new number starting at
620:     //     the top and going downward. I know these might not match up with
621:     //     the iterator order, but we can fix it later.
622:     template<typename Mesh>
623:     static void buildDualCSR(const Obj<Mesh>& mesh, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
624:       const Obj<typename Mesh::sieve_type>&         sieve        = mesh->getSieve();
625:       const Obj<typename Mesh::label_sequence>&     cells        = mesh->heightStratum(0);
626:       const typename Mesh::label_sequence::iterator cEnd         = cells->end();
627:       const int                                     numCells     = cells->size();
628:       int                                           newCell      = numCells;
629:       Obj<typename Mesh::sieve_type>                overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
630:       int                                           offset       = 0;
631:       const int                                     cellOffset   = zeroBase ? 0 : 1;
632:       const int                                     dim          = mesh->getDimension();
633:       std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;

635:       // TODO: This is necessary for parallel partitioning
636:       //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
637:       if (numCells == 0) {
638:         *numVertices = 0;
639:         *offsets     = NULL;
640:         *adjacency   = NULL;
641:         return;
642:       }
643:       int *off = alloc_type().allocate(numCells+1);
644:       int *adj;
645:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
646:       if (mesh->depth() == dim) {
647:         int c = 1;

649:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
650:           const Obj<typename Mesh::sieve_type::traits::coneSequence>&     faces = sieve->cone(*c_iter);
651:           const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd  = faces->end();

653:           off[c] = off[c-1];
654:           for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
655:             if (sieve->support(*f_iter)->size() == 2) {
656:               off[c]++;
657:             } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
658:               off[c]++;
659:             }
660:           }
661:         }
662:         adj = alloc_type().allocate(off[numCells]);
663:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
664:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
665:           const Obj<typename Mesh::sieve_type::traits::coneSequence>&     faces = sieve->cone(*c_iter);
666:           const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd  = faces->end();

668:           for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
669:             const Obj<typename Mesh::sieve_type::traits::supportSequence>&     neighbors = sieve->support(*f_iter);
670:             const typename Mesh::sieve_type::traits::supportSequence::iterator nEnd      = neighbors->end();

672:             for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
673:               if (*n_iter != *c_iter) adj[offset++] = *n_iter + cellOffset;
674:             }
675:             const Obj<typename Mesh::sieve_type::traits::supportSequence>&     oNeighbors = overlapSieve->support(*f_iter);
676:             const typename Mesh::sieve_type::traits::supportSequence::iterator onEnd      = oNeighbors->end();

678:             for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = oNeighbors->begin(); n_iter != onEnd; ++n_iter) {
679:               adj[offset++] = *n_iter + cellOffset;
680:             }
681:           }
682:         }
683:       } else if (mesh->depth() == 1) {
684:         std::set<typename Mesh::point_type> *neighborCells = new std::set<typename Mesh::point_type>[numCells];
685:         const int                            corners       = sieve->cone(*cells->begin())->size();
686:         int                                  faceVertices;

688:         if (corners == dim+1) {
689:           faceVertices = dim;
690:         } else if ((dim == 2) && (corners == 4)) {
691:           faceVertices = 2;
692:         } else if ((dim == 3) && (corners == 8)) {
693:           faceVertices = 4;
694:         } else {
695:           throw ALE::Exception("Could not determine number of face vertices");
696:         }
697:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
698:             const Obj<typename Mesh::sieve_type::traits::coneSequence>&     vertices = sieve->cone(*c_iter);
699:             const typename Mesh::sieve_type::traits::coneSequence::iterator vEnd     = vertices->end();

701:             for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
702:               const Obj<typename Mesh::sieve_type::traits::supportSequence>&     neighbors = sieve->support(*v_iter);
703:               const typename Mesh::sieve_type::traits::supportSequence::iterator nEnd      = neighbors->end();

705:               for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
706:                 if (*c_iter == *n_iter) continue;
707:                 if ((int) sieve->nMeet(*c_iter, *n_iter, 1)->size() == faceVertices) {
708:                   if ((*c_iter < numCells) && (*n_iter < numCells)) {
709:                     neighborCells[*c_iter].insert(*n_iter);
710:                   } else {
711:                     typename Mesh::point_type e = *c_iter, n = *n_iter;

713:                     if (*c_iter >= numCells) {
714:                       if (newCells.find(*c_iter) == newCells.end()) newCells[*c_iter] = --newCell;
715:                       e = newCells[*c_iter];
716:                     }
717:                     if (*n_iter >= numCells) {
718:                       if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
719:                       n = newCells[*n_iter];
720:                     }
721:                     neighborCells[e].insert(n);
722:                   }
723:                 }
724:               }
725:             }
726:           }
727:           off[0] = 0;
728:           for(int c = 1; c <= numCells; c++) {
729:             off[c] = neighborCells[c-1].size() + off[c-1];
730:           }
731:           adj = alloc_type().allocate(off[numCells]);
732:           for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
733:           for(int c = 0; c < numCells; c++) {
734:             for(typename std::set<typename Mesh::point_type>::iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
735:               adj[offset++] = *n_iter + cellOffset;
736:             }
737:           }
738:           delete [] neighborCells;
739:       } else {
740:         throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
741:       }
742:       if (offset != off[numCells]) {
743:         ostringstream msg;
744:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
745:         throw ALE::Exception(msg.str().c_str());
746:       }
747:       *numVertices = numCells;
748:       *offsets     = off;
749:       *adjacency   = adj;
750:     };
751:     template<typename Mesh>
752:     static void buildDualCSRV(const Obj<Mesh>& mesh, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
753:       const Obj<typename Mesh::sieve_type>&         sieve        = mesh->getSieve();
754:       const Obj<typename Mesh::label_sequence>&     cells        = mesh->heightStratum(0);
755:       const typename Mesh::label_sequence::iterator cEnd         = cells->end();
756:       const int                                     numCells     = cells->size();
757:       Obj<typename Mesh::sieve_type>                overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
758:       int                                           offset       = 0;
759:       const int                                     cellOffset   = zeroBase ? 0 : 1;
760:       const int                                     dim          = mesh->getDimension();
761:       std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;

763:       // TODO: This is necessary for parallel partitioning
764:       //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
765:       overlapSieve->setChart(sieve->getChart());
766:       overlapSieve->allocate();
767:       if (numCells == 0) {
768:         *numVertices = 0;
769:         *offsets     = NULL;
770:         *adjacency   = NULL;
771:         return;
772:       }
773:       int *off = alloc_type().allocate(numCells+1);
774:       int *adj;
775:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
776:       if (mesh->depth() == dim) {
777:         OffsetVisitor<typename Mesh::sieve_type> oV(*sieve, *overlapSieve, off);

779:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
780:           sieve->cone(*c_iter, oV);
781:         }
782:         for(int p = 1; p <= numCells; ++p) {
783:           off[p] = off[p] + off[p-1];
784:         }
785:         adj = alloc_type().allocate(off[numCells]);
786:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
787:         AdjVisitor<typename Mesh::sieve_type> aV(adj, zeroBase);
788:         ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > sV(*sieve, aV);
789:         ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > ovSV(*overlapSieve, aV);

791:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
792:           aV.setCell(*c_iter);
793:           sieve->cone(*c_iter, sV);
794:           sieve->cone(*c_iter, ovSV);
795:         }
796:         offset = aV.getOffset();
797:       } else if (mesh->depth() == 1) {
798:         typedef MeetVisitor<typename Mesh::sieve_type>    mv_type;
799:         typedef typename ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, mv_type> sv_type;
800:         const int corners = sieve->getConeSize(*cells->begin());
801:         int       faceVertices;

803:         if (corners == dim+1) {
804:           faceVertices = dim;
805:         } else if ((dim == 2) && (corners == 4)) {
806:           faceVertices = 2;
807:         } else if ((dim == 3) && (corners == 8)) {
808:           faceVertices = 4;
809:         } else {
810:           throw ALE::Exception("Could not determine number of face vertices");
811:         }
812:         mv_type mV(*sieve, numCells, faceVertices);
813:         sv_type sV(*sieve, mV);

815:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
816:           mV.setCell(*c_iter);
817:           sieve->cone(*c_iter, sV);
818:         }
819:         const typename mv_type::neighbors_type *neighborCells = mV.getNeighbors();

821:         off[0] = 0;
822:         for(int c = 1; c <= numCells; c++) {
823:           off[c] = neighborCells[c-1].size() + off[c-1];
824:         }
825:         adj = alloc_type().allocate(off[numCells]);
826:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
827:         for(int c = 0; c < numCells; c++) {
828:           for(typename mv_type::neighbors_type::const_iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
829:             //std::cout << "Adding dual edge from " << c << " to " << *n_iter << std::endl;
830:             adj[offset++] = *n_iter + cellOffset;
831:           }
832:         }
833:       } else {
834:         throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
835:       }
836:       if (offset != off[numCells]) {
837:         ostringstream msg;
838:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
839:         throw ALE::Exception(msg.str().c_str());
840:       }
841:       *numVertices = numCells;
842:       *offsets     = off;
843:       *adjacency   = adj;
844:     };
845:     // This produces a hypergraph (each face is a vertex and each cell is a hyperedge)
846:     //   numbering: A contiguous numbering of the faces
847:     //   numEdges:  The number of edges in the hypergraph
848:     //   adjacency: The vertices in each edge
849:     template<typename Mesh>
850:     static void buildFaceDualCSR(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
851:       const Obj<typename Mesh::sieve_type>&         sieve = mesh->getSieve();
852:       const Obj<typename Mesh::label_sequence>&     cells = mesh->heightStratum(0);
853:       const typename Mesh::label_sequence::iterator cEnd  = cells->end();
854:       const int faceOffset = zeroBase ? 0 : 1;
855:       int       numCells = cells->size();
856:       int       c        = 1;

858:       if (mesh->depth() != mesh->getDimension()) {throw ALE::Exception("Not yet implemented for non-interpolated meshes");}
859:       int *off = alloc_type().allocate(numCells+1);
860:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
861:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
862:         off[c] = sieve->cone(*c_iter)->size() + off[c-1];
863:       }
864:       int *adj = alloc_type().allocate(off[numCells]);
865:       for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
866:       int  offset = 0;
867:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
868:         const Obj<typename Mesh::sieve_type::traits::coneSequence>&     faces = sieve->cone(*c_iter);
869:         const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd  = faces->end();

871:         for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter, ++offset) {
872:           adj[offset] = numbering->getIndex(*f_iter) + faceOffset;
873:         }
874:       }
875:       if (offset != off[numCells]) {
876:         ostringstream msg;
877:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
878:         throw ALE::Exception(msg.str().c_str());
879:       }
880:       *numEdges  = numCells;
881:       *offsets   = off;
882:       *adjacency = adj;
883:     };
884:     template<typename Mesh>
885:     static void buildFaceDualCSRV(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
886:       throw ALE::Exception("Not implemented");
887:     };
888:     static void destroyCSR(int numPoints, int *offsets, int *adjacency) {
889:       if (adjacency) {
890:         for(int i = 0; i < offsets[numPoints]; ++i) {alloc_type().destroy(adjacency+i);}
891:         alloc_type().deallocate(adjacency, offsets[numPoints]);
892:       }
893:       if (offsets) {
894:         for(int i = 0; i < numPoints+1; ++i) {alloc_type().destroy(offsets+i);}
895:         alloc_type().deallocate(offsets, numPoints+1);
896:       }
897:     };
898:     template<typename OldSection, typename Partition, typename Renumbering, typename NewSection>
899:     static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<NewSection>& newSection) {
900:       const typename Partition::value_type *points    = partition->restrictPoint(oldSection->commRank());
901:       const int                             numPoints = partition->getFiberDimension(oldSection->commRank());

903:       for(int p = 0; p < numPoints; ++p) {
904:         if (oldSection->hasPoint(points[p])) {
905:           newSection->setFiberDimension(renumbering[points[p]], oldSection->getFiberDimension(points[p]));
906:         }
907:       }
908:       if (numPoints) {newSection->allocatePoint();}
909:       for(int p = 0; p < numPoints; ++p) {
910:         if (oldSection->hasPoint(points[p])) {
911:           newSection->updatePointAll(renumbering[points[p]], oldSection->restrictPoint(points[p]));
912:         }
913:       }
914:     };
915:     // Specialize to ArrowSection
916:     template<typename OldSection, typename Partition, typename Renumbering>
917:     static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<UniformSection<MinimalArrow<int,int>,int> >& newSection) {
918:       typedef UniformSection<MinimalArrow<int,int>,int> NewSection;
919:       const typename Partition::value_type    *points    = partition->restrictPoint(oldSection->commRank());
920:       const int                                numPoints = partition->getFiberDimension(oldSection->commRank());
921:       const typename OldSection::chart_type&   oldChart  = oldSection->getChart();
922:       std::set<typename Partition::value_type> myPoints;

924:       for(int p = 0; p < numPoints; ++p) {
925:         myPoints.insert(points[p]);
926:       }
927:       for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
928:         if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
929:           newSection->setFiberDimension(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), oldSection->getFiberDimension(*c_iter));
930:         }
931:       }
932:       if (oldChart.size()) {newSection->allocatePoint();}
933:       for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
934:         if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
935:           const typename OldSection::value_type *values = oldSection->restrictPoint(*c_iter);

937:           newSection->updatePointAll(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), values);
938:         }
939:       }
940:     };
941:     template<typename Sifter, typename Section, typename Renumbering>
942:     static void createLocalSifter(const Obj<Sifter>& sifter, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sifter>& localSifter) {
943:       const typename Section::value_type *points    = partition->restrictPoint(sifter->commRank());
944:       const int                           numPoints = partition->getFiberDimension(sifter->commRank());

946:       for(int p = 0; p < numPoints; ++p) {
947:         const Obj<typename Sifter::traits::coneSequence>&     cone = sifter->cone(points[p]);
948:         const typename Sifter::traits::coneSequence::iterator cEnd = cone->end();

950:         for(typename Sifter::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
951:           localSifter->addArrow(*c_iter, renumbering[points[p]]);
952:         }
953:       }
954:     };
955:     template<typename Sieve, typename Section, typename Renumbering>
956:     static void createLocalSieve(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
957:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
958:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());

960:       for(int p = 0; p < numPoints; ++p) {
961:         Obj<typename Sieve::coneSet> current = new typename Sieve::coneSet();
962:         Obj<typename Sieve::coneSet> next    = new typename Sieve::coneSet();
963:         Obj<typename Sieve::coneSet> tmp;

965:         current->insert(points[p]);
966:         while(current->size()) {
967:           for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
968:             const Obj<typename Sieve::traits::coneSequence>& cone = sieve->cone(*p_iter);
969: 
970:             for(typename Sieve::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
971:               localSieve->addArrow(renumbering[*c_iter], renumbering[*p_iter], c_iter.color());
972:               next->insert(*c_iter);
973:             }
974:           }
975:           tmp = current; current = next; next = tmp;
976:           next->clear();
977:         }
978:         if (height) {
979:           current->insert(points[p]);
980:           while(current->size()) {
981:             for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
982:               const Obj<typename Sieve::traits::supportSequence>& support = sieve->support(*p_iter);
983: 
984:               for(typename Sieve::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
985:                 localSieve->addArrow(renumbering[*p_iter], renumbering[*s_iter], s_iter.color());
986:                 next->insert(*s_iter);
987:               }
988:             }
989:             tmp = current; current = next; next = tmp;
990:             next->clear();
991:           }
992:         }
993:       }
994:     };
995:     template<typename Mesh, typename Section, typename Renumbering>
996:     static void createLocalMesh(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
997:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
998:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1000:       createLocalSieve(sieve, partition, renumbering, localSieve, height);
1001:     };
1002:     template<typename Sieve, typename Section, typename Renumbering>
1003:     static void sizeLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1004:       typedef std::set<typename Sieve::point_type> pointSet;
1005:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1006:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1007:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1008:       const pointSet                      pSet(points, &points[numPoints]);
1009:       ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> fV(pSet, renumbering, maxSize);

1011:       for(int p = 0; p < numPoints; ++p) {
1012:         sieve->cone(points[p], fV);
1013:         localSieve->setConeSize(renumbering[points[p]], fV.getSize());
1014:         fV.clear();
1015:         sieve->support(points[p], fV);
1016:         localSieve->setSupportSize(renumbering[points[p]], fV.getSize());
1017:         fV.clear();
1018:       }
1019:     };
1020:     template<typename Mesh, typename Section, typename Renumbering>
1021:     static void sizeLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1022:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
1023:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1025:       sizeLocalSieveV(sieve, partition, renumbering, localSieve, height);
1026:     };
1027:     template<typename Sieve, typename Section, typename Renumbering>
1028:     static void createLocalLabelV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1029:       typedef std::set<typename Sieve::point_type> pointSet;
1030:       typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1031:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1032:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1033:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1034:       typename Sieve::point_type         *oPoints   = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1035:       int                                *oOrients  = new int[std::max(1, sieve->getMaxConeSize())];
1036:       const pointSet                      pSet(points, &points[numPoints]);
1037:       visitor_type fV(pSet, renumbering, maxSize);

1039:       for(int p = 0; p < numPoints; ++p) {
1040:         fV.useRenumbering(false);
1041:         sieve->orientedCone(points[p], fV);
1042:         const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1043:         const int                                         n = fV.getOrientedSize();
1044:         for(int i = 0; i < n; ++i) {
1045:           oPoints[i]  = q[i].first;
1046:           oOrients[i] = q[i].second;
1047:         }
1048:         localSieve->setCone(oPoints, renumbering[points[p]]);
1049:         localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1050:         fV.clear();
1051:         fV.useRenumbering(true);
1052:         sieve->support(points[p], fV);
1053:         if (fV.getSize()) {localSieve->setSupport(points[p], fV.getPoints());}
1054:         fV.clear();
1055:       }
1056:       delete [] oPoints;
1057:       delete [] oOrients;
1058:     };
1059:     template<typename Sieve, typename Section, typename Renumbering>
1060:     static void createLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1061:       typedef std::set<typename Sieve::point_type> pointSet;
1062:       typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1063:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1064:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1065:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1066:       typename Sieve::point_type         *oPoints   = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1067:       int                                *oOrients  = new int[std::max(1, sieve->getMaxConeSize())];
1068:       const pointSet                      pSet(points, &points[numPoints]);
1069:       visitor_type fV(pSet, renumbering, maxSize);

1071:       for(int p = 0; p < numPoints; ++p) {
1072:         ///sieve->cone(points[p], fV);
1073:         ///localSiaeve->setCone(fV.getPoints(), renumbering[points[p]]);
1074:         sieve->orientedCone(points[p], fV);
1075:         const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1076:         const int                                         n = fV.getOrientedSize();
1077:         for(int i = 0; i < n; ++i) {
1078:           oPoints[i]  = q[i].first;
1079:           oOrients[i] = q[i].second;
1080:         }
1081:         localSieve->setCone(oPoints, renumbering[points[p]]);
1082:         localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1083:         fV.clear();
1084:         sieve->support(points[p], fV);
1085:         localSieve->setSupport(renumbering[points[p]], fV.getPoints());
1086:         fV.clear();
1087:       }
1088:       delete [] oPoints;
1089:       delete [] oOrients;
1090:     };
1091:     template<typename Mesh, typename Section, typename Renumbering>
1092:     static void createLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1093:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
1094:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1096:       createLocalSieveV(sieve, partition, renumbering, localSieve, height);
1097:     };
1098:   public: // Partitioning
1099:     //   partition:    Should be properly allocated on input
1100:     //   height:       Height of the point set to uniquely partition
1101:     // TODO: Could rebind assignment section to the type of the output
1102:     template<typename Mesh, typename Section>
1103:     static void createPartition(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1104:       MeshManager<Mesh> manager(mesh, height);
1105:       int              *start     = NULL;
1106:       int              *adjacency = NULL;

1108:       if (height == 0) {
1109:         int numVertices;

1111:         buildDualCSR(mesh, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1112:         GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1113:         destroyCSR(numVertices, start, adjacency);
1114:       } else if (height == 1) {
1115:         int numEdges;

1117:         buildFaceDualCSR(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1118:         GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1119:         destroyCSR(numEdges, start, adjacency);
1120:       } else {
1121:         throw ALE::Exception("Invalid partition height");
1122:       }
1123:     };
1124:     template<typename Mesh, typename Section>
1125:     static void createPartitionV(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1126:       MeshManager<Mesh> manager(mesh, height);
1127:       int              *start     = NULL;
1128:       int              *adjacency = NULL;

1130:       PETSc::Log::Event("PartitionCreate").begin();
1131:       if (height == 0) {
1132:         int numVertices;

1134:         buildDualCSRV(mesh, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1135:         GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1136:         destroyCSR(numVertices, start, adjacency);
1137:       } else if (height == 1) {
1138:         int numEdges;

1140:         throw ALE::Exception("Not yet implemented");
1141: #if 0
1142:         buildFaceDualCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1143: #endif
1144:         GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1145:         destroyCSR(numEdges, start, adjacency);
1146:       } else {
1147:         throw ALE::Exception("Invalid partition height");
1148:       }
1149:       PETSc::Log::Event("PartitionCreate").end();
1150:     };
1151:     // Add in the points in the closure (and star) of the partitioned points
1152:     template<typename Mesh, typename Section>
1153:     static void createPartitionClosure(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1154:       const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1155:       const typename Section::chart_type&   chart = pointPartition->getChart();
1156:       size_t                                size  = 0;

1158:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1159:         const typename Section::value_type    *points    = pointPartition->restrictPoint(*r_iter);
1160:         const int                              numPoints = pointPartition->getFiberDimension(*r_iter);
1161:         std::set<typename Section::value_type> closure;

1163:         // TODO: Use Quiver's closure() here instead
1164:         for(int p = 0; p < numPoints; ++p) {
1165:           Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1166:           Obj<typename Mesh::sieve_type::coneSet> next    = new typename Mesh::sieve_type::coneSet();
1167:           Obj<typename Mesh::sieve_type::coneSet> tmp;

1169:           current->insert(points[p]);
1170:           closure.insert(points[p]);
1171:           while(current->size()) {
1172:             for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1173:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1174: 
1175:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1176:                 closure.insert(*c_iter);
1177:                 next->insert(*c_iter);
1178:               }
1179:             }
1180:             tmp = current; current = next; next = tmp;
1181:             next->clear();
1182:           }
1183:           if (height) {
1184:             current->insert(points[p]);
1185:             while(current->size()) {
1186:               for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1187:                 const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1188: 
1189:                 for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1190:                   closure.insert(*s_iter);
1191:                   next->insert(*s_iter);
1192:                 }
1193:               }
1194:               tmp = current; current = next; next = tmp;
1195:               next->clear();
1196:             }
1197:           }
1198:         }
1199:         partition->setFiberDimension(*r_iter, closure.size());
1200:         size = std::max(size, closure.size());
1201:       }
1202:       partition->allocatePoint();
1203:       typename Section::value_type *values = new typename Section::value_type[size];

1205:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1206:         const typename Section::value_type    *points    = pointPartition->restrictPoint(*r_iter);
1207:         const int                              numPoints = pointPartition->getFiberDimension(*r_iter);
1208:         std::set<typename Section::value_type> closure;

1210:         // TODO: Use Quiver's closure() here instead
1211:         for(int p = 0; p < numPoints; ++p) {
1212:           Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1213:           Obj<typename Mesh::sieve_type::coneSet> next    = new typename Mesh::sieve_type::coneSet();
1214:           Obj<typename Mesh::sieve_type::coneSet> tmp;

1216:           current->insert(points[p]);
1217:           closure.insert(points[p]);
1218:           while(current->size()) {
1219:             for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1220:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1221: 
1222:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1223:                 closure.insert(*c_iter);
1224:                 next->insert(*c_iter);
1225:               }
1226:             }
1227:             tmp = current; current = next; next = tmp;
1228:             next->clear();
1229:           }
1230:           if (height) {
1231:             current->insert(points[p]);
1232:             while(current->size()) {
1233:               for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1234:                 const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1235: 
1236:                 for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1237:                   closure.insert(*s_iter);
1238:                   next->insert(*s_iter);
1239:                 }
1240:               }
1241:               tmp = current; current = next; next = tmp;
1242:               next->clear();
1243:             }
1244:           }
1245:         }
1246:         int i = 0;

1248:         for(typename std::set<typename Section::value_type>::const_iterator p_iter = closure.begin(); p_iter != closure.end(); ++p_iter, ++i) {
1249:           values[i] = *p_iter;
1250:         }
1251:         partition->updatePoint(*r_iter, values);
1252:       }
1253:       delete [] values;
1254:     };
1255:     template<typename Mesh, typename Section>
1256:     static void createPartitionClosureV(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1257:       typedef ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type> visitor_type;
1258:       const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1259:       const typename Section::chart_type&   chart = pointPartition->getChart();
1260:       size_t                                size  = 0;

1262:       PETSc::Log::Event("PartitionClosure").begin();
1263:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1264:         const typename Section::value_type *points    = pointPartition->restrictPoint(*r_iter);
1265:         const int                           numPoints = pointPartition->getFiberDimension(*r_iter);
1266:         typename visitor_type::visitor_type nV;
1267:         visitor_type                        cV(*sieve, nV);

1269:         for(int p = 0; p < numPoints; ++p) {
1270:           sieve->cone(points[p], cV);
1271:           if (height) {
1272:             cV.setIsCone(false);
1273:             sieve->support(points[p], cV);
1274:           }
1275:         }
1276:         partition->setFiberDimension(*r_iter, cV.getPoints().size());
1277:         size = std::max(size, cV.getPoints().size());
1278:       }
1279:       partition->allocatePoint();
1280:       typename Section::value_type *values = new typename Section::value_type[size];

1282:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1283:         const typename Section::value_type *points    = pointPartition->restrictPoint(*r_iter);
1284:         const int                           numPoints = pointPartition->getFiberDimension(*r_iter);
1285:         typename visitor_type::visitor_type nV;
1286:         visitor_type                        cV(*sieve, nV);

1288:         for(int p = 0; p < numPoints; ++p) {
1289:           sieve->cone(points[p], cV);
1290:           if (height) {
1291:             cV.setIsCone(false);
1292:             sieve->support(points[p], cV);
1293:           }
1294:         }
1295:         int i = 0;

1297:         for(typename std::set<typename Mesh::point_type>::const_iterator p_iter = cV.getPoints().begin(); p_iter != cV.getPoints().end(); ++p_iter, ++i) {
1298:           values[i] = *p_iter;
1299:         }
1300:         partition->updatePoint(*r_iter, values);
1301:       }
1302:       delete [] values;
1303:       PETSc::Log::Event("PartitionClosure").end();
1304:     };
1305:     // Create a section mapping points to partitions
1306:     template<typename Section, typename MapSection>
1307:     static void createPartitionMap(const Obj<Section>& partition, const Obj<MapSection>& partitionMap) {
1308:       const typename Section::chart_type& chart = partition->getChart();

1310:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1311:         partitionMap->setFiberDimension(*p_iter, 1);
1312:       }
1313:       partitionMap->allocatePoint();
1314:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1315:         const typename Section::value_type *points = partition->restrictPoint(*p_iter);
1316:         const int                           size   = partition->getFiberDimension(*p_iter);
1317:         const typename Section::point_type  part   = *p_iter;

1319:         for(int i = 0; i < size; ++i) {
1320:           partitionMap->updatePoint(points[i], &part);
1321:         }
1322:       }
1323:     };
1324:   };
1325: #endif

1327:   namespace New {
1328:     template<typename Bundle_, typename Alloc_ = typename Bundle_::alloc_type>
1329:     class Partitioner {
1330:     public:
1331:       typedef Bundle_                          bundle_type;
1332:       typedef Alloc_                           alloc_type;
1333:       typedef typename bundle_type::sieve_type sieve_type;
1334:       typedef typename bundle_type::point_type point_type;
1335:     public:
1338:       // This creates a CSR representation of the adjacency matrix for cells
1339:       // - We allow an exception to contiguous numbering.
1340:       //   If the cell id > numElements, we assign a new number starting at
1341:       //     the top and going downward. I know these might not match up with
1342:       //     the iterator order, but we can fix it later.
1343:       static void buildDualCSR(const Obj<bundle_type>& bundle, const int dim, int **offsets, int **adjacency) {
1344:         ALE_LOG_EVENT_BEGIN;
1345:         typedef typename ALE::New::Completion<bundle_type, point_type, alloc_type> completion;
1346:         const Obj<sieve_type>&                           sieve        = bundle->getSieve();
1347:         const Obj<typename bundle_type::label_sequence>& elements     = bundle->heightStratum(0);
1348:         Obj<sieve_type>                                  overlapSieve = new sieve_type(bundle->comm(), bundle->debug());
1349:         std::map<point_type, point_type>                 newCells;
1350:         int  numElements = elements->size();
1351:         int  newCell     = numElements;
1352:         int *off         = new int[numElements+1];
1353:         int  offset      = 0;
1354:         int *adj;

1356:         completion::scatterSupports(sieve, overlapSieve, bundle->getSendOverlap(), bundle->getRecvOverlap(), bundle);
1357:         if (numElements == 0) {
1358:           *offsets   = NULL;
1359:           *adjacency = NULL;
1360:           ALE_LOG_EVENT_END;
1361:           return;
1362:         }
1363:         if (bundle->depth() == dim) {
1364:           int e = 1;

1366:           off[0] = 0;
1367:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1368:             const Obj<typename sieve_type::traits::coneSequence>& faces  = sieve->cone(*e_iter);
1369:             typename sieve_type::traits::coneSequence::iterator   fBegin = faces->begin();
1370:             typename sieve_type::traits::coneSequence::iterator   fEnd   = faces->end();

1372:             off[e] = off[e-1];
1373:             for(typename sieve_type::traits::coneSequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1374:               if (sieve->support(*f_iter)->size() == 2) {
1375:                 off[e]++;
1376:               } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
1377:                 off[e]++;
1378:               }
1379:             }
1380:             e++;
1381:           }
1382:           adj = new int[off[numElements]];
1383:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1384:             const Obj<typename sieve_type::traits::coneSequence>& faces  = sieve->cone(*e_iter);
1385:             typename sieve_type::traits::coneSequence::iterator   fBegin = faces->begin();
1386:             typename sieve_type::traits::coneSequence::iterator   fEnd   = faces->end();

1388:             for(typename sieve_type::traits::coneSequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1389:               const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*f_iter);
1390:               typename sieve_type::traits::supportSequence::iterator   nBegin    = neighbors->begin();
1391:               typename sieve_type::traits::supportSequence::iterator   nEnd      = neighbors->end();

1393:               for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
1394:                 if (*n_iter != *e_iter) adj[offset++] = *n_iter;
1395:               }
1396:               const Obj<typename sieve_type::traits::supportSequence>& oNeighbors = overlapSieve->support(*f_iter);
1397:               typename sieve_type::traits::supportSequence::iterator   onBegin    = oNeighbors->begin();
1398:               typename sieve_type::traits::supportSequence::iterator   onEnd      = oNeighbors->end();

1400:               for(typename sieve_type::traits::supportSequence::iterator n_iter = onBegin; n_iter != onEnd; ++n_iter) {
1401:                 adj[offset++] = *n_iter;
1402:               }
1403:             }
1404:           }
1405:         } else if (bundle->depth() == 1) {
1406:           std::set<point_type> *neighborCells = new std::set<point_type>[numElements];
1407:           int corners      = sieve->cone(*elements->begin())->size();
1408:           int faceVertices = -1;

1410:           if (corners == dim+1) {
1411:             faceVertices = dim;
1412:           } else if ((dim == 2) && (corners == 4)) {
1413:             faceVertices = 2;
1414:           } else if ((dim == 3) && (corners == 8)) {
1415:             faceVertices = 4;
1416:           } else {
1417:             throw ALE::Exception("Could not determine number of face vertices");
1418:           }
1419:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1420:             const Obj<typename sieve_type::traits::coneSequence>& vertices  = sieve->cone(*e_iter);
1421:             typename sieve_type::traits::coneSequence::iterator vEnd = vertices->end();

1423:             for(typename sieve_type::traits::coneSequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1424:               const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
1425:               typename sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();

1427:               for(typename sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
1428:                 if (*e_iter == *n_iter) continue;
1429:                 if ((int) sieve->nMeet(*e_iter, *n_iter, 1)->size() == faceVertices) {
1430:                   if ((*e_iter < numElements) && (*n_iter < numElements)) {
1431:                     neighborCells[*e_iter].insert(*n_iter);
1432:                   } else {
1433:                     point_type e = *e_iter, n = *n_iter;

1435:                     if (*e_iter >= numElements) {
1436:                       if (newCells.find(*e_iter) == newCells.end()) newCells[*e_iter] = --newCell;
1437:                       e = newCells[*e_iter];
1438:                     }
1439:                     if (*n_iter >= numElements) {
1440:                       if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
1441:                       n = newCells[*n_iter];
1442:                     }
1443:                     neighborCells[e].insert(n);
1444:                   }
1445:                 }
1446:               }
1447:             }
1448:           }
1449:           off[0] = 0;
1450:           for(int e = 1; e <= numElements; e++) {
1451:             off[e] = neighborCells[e-1].size() + off[e-1];
1452:           }
1453:           adj = new int[off[numElements]];
1454:           for(int e = 0; e < numElements; e++) {
1455:             for(typename std::set<point_type>::iterator n_iter = neighborCells[e].begin(); n_iter != neighborCells[e].end(); ++n_iter) {
1456:               adj[offset++] = *n_iter;
1457:             }
1458:           }
1459:           delete [] neighborCells;
1460:         } else {
1461:           throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
1462:         }
1463:         if (offset != off[numElements]) {
1464:           ostringstream msg;
1465:           msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1466:           throw ALE::Exception(msg.str().c_str());
1467:         }
1468:         //std::cout << "numElements: " << numElements << " newCell: " << newCell << std::endl;
1469:         *offsets   = off;
1470:         *adjacency = adj;
1471:         ALE_LOG_EVENT_END;
1472:       };
1475:       // This creates a CSR representation of the adjacency hypergraph for faces
1476:       static void buildFaceCSR(const Obj<bundle_type>& bundle, const int dim, const Obj<typename bundle_type::numbering_type>& fNumbering, int *numEdges, int **offsets, int **adjacency) {
1477:         ALE_LOG_EVENT_BEGIN;
1478:         const Obj<sieve_type>&                           sieve    = bundle->getSieve();
1479:         const Obj<typename bundle_type::label_sequence>& elements = bundle->heightStratum(0);
1480:         int  numElements = elements->size();
1481:         int *off         = new int[numElements+1];
1482:         int  e;

1484:         if (bundle->depth() != dim) {
1485:           throw ALE::Exception("Not yet implemented for non-interpolated meshes");
1486:         }
1487:         off[0] = 0;
1488:         e      = 1;
1489:         for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1490:           off[e] = sieve->cone(*e_iter)->size() + off[e-1];
1491:           e++;
1492:         }
1493:         int *adj    = new int[off[numElements]];
1494:         int  offset = 0;
1495:         for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1496:           const Obj<typename sieve_type::traits::coneSequence>& faces = sieve->cone(*e_iter);
1497:           typename sieve_type::traits::coneSequence::iterator   fEnd  = faces->end();

1499:           for(typename sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
1500:             adj[offset++] = fNumbering->getIndex(*f_iter);
1501:           }
1502:         }
1503:         if (offset != off[numElements]) {
1504:           ostringstream msg;
1505:           msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1506:           throw ALE::Exception(msg.str().c_str());
1507:         }
1508:         *numEdges  = numElements;
1509:         *offsets   = off;
1510:         *adjacency = adj;
1511:         ALE_LOG_EVENT_END;
1512:       };
1513:       template<typename PartitionType>
1514:       static PartitionType *subordinatePartition(const Obj<bundle_type>& bundle, int levels, const Obj<bundle_type>& subBundle, const PartitionType assignment[]) {
1515:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth());
1516:         const Obj<typename bundle_type::label_sequence>& cells      = subBundle->heightStratum(0);
1517:         const Obj<typename bundle_type::numbering_type>& sNumbering = bundle->getFactory()->getLocalNumbering(subBundle, subBundle->depth());
1518:         const int        numCells      = cells->size();
1519:         PartitionType   *subAssignment = new PartitionType[numCells];

1521:         if (levels != 1) {
1522:           throw ALE::Exception("Cannot calculate subordinate partition for any level separation other than 1");
1523:         } else {
1524:           const Obj<typename bundle_type::sieve_type>&   sieve    = bundle->getSieve();
1525:           const Obj<typename bundle_type::sieve_type>&   subSieve = subBundle->getSieve();
1526:           Obj<typename bundle_type::sieve_type::coneSet> tmpSet   = new typename bundle_type::sieve_type::coneSet();
1527:           Obj<typename bundle_type::sieve_type::coneSet> tmpSet2  = new typename bundle_type::sieve_type::coneSet();

1529:           for(typename bundle_type::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1530:             const Obj<typename bundle_type::sieve_type::coneSequence>& cone = subSieve->cone(*c_iter);

1532:             Obj<typename bundle_type::sieve_type::supportSet> cell = sieve->nJoin1(cone);
1533:             if (cell->size() != 1) {
1534:               std::cout << "Indeterminate subordinate partition for face " << *c_iter << std::endl;
1535:               for(typename bundle_type::sieve_type::supportSet::iterator s_iter = cell->begin(); s_iter != cell->end(); ++s_iter) {
1536:                 std::cout << "  cell " << *s_iter << std::endl;
1537:               }
1538:               // Could relax this to choosing the first one
1539:               throw ALE::Exception("Indeterminate subordinate partition");
1540:             }
1541:             subAssignment[sNumbering->getIndex(*c_iter)] = assignment[cNumbering->getIndex(*cell->begin())];
1542:             tmpSet->clear();
1543:             tmpSet2->clear();
1544:           }
1545:         }
1546:         return subAssignment;
1547:       };
1548:     };
1549: #ifdef PETSC_HAVE_CHACO
1550:     namespace Chaco {
1551:       template<typename Bundle_>
1552:       class Partitioner {
1553:       public:
1554:         typedef Bundle_                          bundle_type;
1555:         typedef typename bundle_type::sieve_type sieve_type;
1556:         typedef typename bundle_type::point_type point_type;
1557:         typedef short int                        part_type;
1558:       public:
1561:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1562:           part_type *assignment = NULL; /* set number of each vtx (length n) */
1563:           int       *start;             /* start of edge list for each vertex */
1564:           int       *adjacency;         /* = adj -> j; edge list data  */

1566:           ALE_LOG_EVENT_BEGIN;
1567:           ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &start, &adjacency);
1568:           if (bundle->commRank() == 0) {
1569:             /* arguments for Chaco library */
1570:             FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1571:             int nvtxs;                              /* number of vertices in full graph */
1572:             int *vwgts = NULL;                      /* weights for all vertices */
1573:             float *ewgts = NULL;                    /* weights for all edges */
1574:             float *x = NULL, *y = NULL, *z = NULL;  /* coordinates for inertial method */
1575:             char *outassignname = NULL;             /*  name of assignment output file */
1576:             char *outfilename = NULL;               /* output file name */
1577:             int architecture = dim;                 /* 0 => hypercube, d => d-dimensional mesh */
1578:             int ndims_tot = 0;                      /* total number of cube dimensions to divide */
1579:             int mesh_dims[3];                       /* dimensions of mesh of processors */
1580:             double *goal = NULL;                    /* desired set sizes for each set */
1581:             int global_method = 1;                  /* global partitioning algorithm */
1582:             int local_method = 1;                   /* local partitioning algorithm */
1583:             int rqi_flag = 0;                       /* should I use RQI/Symmlq eigensolver? */
1584:             int vmax = 200;                         /* how many vertices to coarsen down to? */
1585:             int ndims = 1;                          /* number of eigenvectors (2^d sets) */
1586:             double eigtol = 0.001;                  /* tolerance on eigenvectors */
1587:             long seed = 123636512;                  /* for random graph mutations */
1588:             float *vCoords[3];

1591:             PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_global_method", &global_method, PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1592:             PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_local_method",  &local_method,  PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1593:             if (global_method == 3) {
1594:               // Inertial Partitioning
1595:               PetscMalloc3(nvtxs,float,&x,nvtxs,float,&y,nvtxs,float,&z);CHKERROR(ierr, "Error in PetscMalloc");
1596:               vCoords[0] = x; vCoords[1] = y; vCoords[2] = z;
1597:               const Obj<typename bundle_type::label_sequence>&    cells       = bundle->heightStratum(0);
1598:               const Obj<typename bundle_type::real_section_type>& coordinates = bundle->getRealSection("coordinates");
1599:               const int corners = bundle->size(coordinates, *(cells->begin()))/dim;
1600:               int       c       = 0;

1602:               for(typename bundle_type::label_sequence::iterator c_iter = cells->begin(); c_iter !=cells->end(); ++c_iter, ++c) {
1603:                 const double *coords = bundle->restrictClosure(coordinates, *c_iter);

1605:                 for(int d = 0; d < dim; ++d) {
1606:                   vCoords[d][c] = 0.0;
1607:                 }
1608:                 for(int v = 0; v < corners; ++v) {
1609:                   for(int d = 0; d < dim; ++d) {
1610:                     vCoords[d][c] += coords[v*dim+d];
1611:                   }
1612:                 }
1613:                 for(int d = 0; d < dim; ++d) {
1614:                   vCoords[d][c] /= corners;
1615:                 }
1616:               }
1617:             }

1619:             nvtxs = bundle->heightStratum(0)->size();
1620:             mesh_dims[0] = bundle->commSize(); mesh_dims[1] = 1; mesh_dims[2] = 1;
1621:             for(int e = 0; e < start[nvtxs]; e++) {
1622:               adjacency[e]++;
1623:             }
1624:             assignment = new part_type[nvtxs];
1625:             PetscMemzero(assignment, nvtxs * sizeof(part_type));CHKERROR(ierr, "Error in PetscMemzero");

1627:             /* redirect output to buffer: chaco -> msgLog */
1628: #ifdef PETSC_HAVE_UNISTD_H
1629:             char *msgLog;
1630:             int fd_stdout, fd_pipe[2], count;

1632:             fd_stdout = dup(1);
1633:             pipe(fd_pipe);
1634:             close(1);
1635:             dup2(fd_pipe[1], 1);
1636:             msgLog = new char[16284];
1637: #endif

1639:             interface(nvtxs, start, adjacency, vwgts, ewgts, x, y, z,
1640:                              outassignname, outfilename, assignment, architecture, ndims_tot,
1641:                              mesh_dims, goal, global_method, local_method, rqi_flag, vmax, ndims,
1642:                              eigtol, seed);

1644: #ifdef PETSC_HAVE_UNISTD_H
1645:             int SIZE_LOG  = 10000;

1647:             fflush(stdout);
1648:             count = read(fd_pipe[0], msgLog, (SIZE_LOG - 1) * sizeof(char));
1649:             if (count < 0) count = 0;
1650:             msgLog[count] = 0;
1651:             close(1);
1652:             dup2(fd_stdout, 1);
1653:             close(fd_stdout);
1654:             close(fd_pipe[0]);
1655:             close(fd_pipe[1]);
1656:             if (bundle->debug()) {
1657:               std::cout << msgLog << std::endl;
1658:             }
1659:             delete [] msgLog;
1660: #endif
1661:             if (global_method == 3) {
1662:               // Inertial Partitioning
1663:               PetscFree3(x, y, z);CHKERROR(ierr, "Error in PetscFree");
1664:             }
1665:           }
1666:           if (adjacency) delete [] adjacency;
1667:           if (start)     delete [] start;
1668:           ALE_LOG_EVENT_END;
1669:           return assignment;
1670:         };
1671:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1672:           throw ALE::Exception("Chaco cannot partition a mesh by faces");
1673:         };
1674:       };
1675:     };
1676: #endif
1677: #ifdef PETSC_HAVE_PARMETIS
1678:     namespace ParMetis {
1679:       template<typename Bundle_>
1680:       class Partitioner {
1681:       public:
1682:         typedef Bundle_                          bundle_type;
1683:         typedef typename bundle_type::sieve_type sieve_type;
1684:         typedef typename bundle_type::point_type point_type;
1685:         typedef int                              part_type;
1686:       public:
1689:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1690:           int    nvtxs      = 0;    // The number of vertices in full graph
1691:           int   *vtxdist;           // Distribution of vertices across processes
1692:           int   *xadj;              // Start of edge list for each vertex
1693:           int   *adjncy;            // Edge lists for all vertices
1694:           int   *vwgt       = NULL; // Vertex weights
1695:           int   *adjwgt     = NULL; // Edge weights
1696:           int    wgtflag    = 0;    // Indicates which weights are present
1697:           int    numflag    = 0;    // Indicates initial offset (0 or 1)
1698:           int    ncon       = 1;    // The number of weights per vertex
1699:           int    nparts     = bundle->commSize(); // The number of partitions
1700:           float *tpwgts;            // The fraction of vertex weights assigned to each partition
1701:           float *ubvec;             // The balance intolerance for vertex weights
1702:           int    options[5];        // Options
1703:           // Outputs
1704:           int    edgeCut;           // The number of edges cut by the partition
1705:           int   *assignment = NULL; // The vertex partition

1707:           options[0] = 0; // Use all defaults
1708:           vtxdist    = new int[nparts+1];
1709:           vtxdist[0] = 0;
1710:           tpwgts     = new float[ncon*nparts];
1711:           for(int p = 0; p < nparts; ++p) {
1712:             tpwgts[p] = 1.0/nparts;
1713:           }
1714:           ubvec      = new float[ncon];
1715:           ubvec[0]   = 1.05;
1716:           nvtxs      = bundle->heightStratum(0)->size();
1717:           assignment = new part_type[nvtxs];
1718:           MPI_Allgather(&nvtxs, 1, MPI_INT, &vtxdist[1], 1, MPI_INT, bundle->comm());
1719:           for(int p = 2; p <= nparts; ++p) {
1720:             vtxdist[p] += vtxdist[p-1];
1721:           }
1722:           if (bundle->commSize() == 1) {
1723:             PetscMemzero(assignment, nvtxs * sizeof(part_type));
1724:           } else {
1725:             ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &xadj, &adjncy);

1727:             if (bundle->debug() && nvtxs) {
1728:               for(int p = 0; p <= nvtxs; ++p) {
1729:                 std::cout << "["<<bundle->commRank()<<"]xadj["<<p<<"] = " << xadj[p] << std::endl;
1730:               }
1731:               for(int i = 0; i < xadj[nvtxs]; ++i) {
1732:                 std::cout << "["<<bundle->commRank()<<"]adjncy["<<i<<"] = " << adjncy[i] << std::endl;
1733:               }
1734:             }
1735:             if (vtxdist[1] == vtxdist[nparts]) {
1736:               if (bundle->commRank() == 0) {
1737:                 METIS_PartGraphKway(&nvtxs, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &nparts, options, &edgeCut, assignment);
1738:                 if (bundle->debug()) {std::cout << "Metis: edgecut is " << edgeCut << std::endl;}
1739:               }
1740:             } else {
1741:               MPI_Comm comm = bundle->comm();

1743:               ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
1744:               if (bundle->debug()) {std::cout << "ParMetis: edgecut is " << edgeCut << std::endl;}
1745:             }
1746:             if (xadj   != NULL) delete [] xadj;
1747:             if (adjncy != NULL) delete [] adjncy;
1748:           }
1749:           delete [] vtxdist;
1750:           delete [] tpwgts;
1751:           delete [] ubvec;
1752:           return assignment;
1753:         };
1756:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1757: #ifdef PETSC_HAVE_HMETIS
1758:           int   *assignment = NULL; // The vertex partition
1759:           int    nvtxs;      // The number of vertices
1760:           int    nhedges;    // The number of hyperedges
1761:           int   *vwgts;      // The vertex weights
1762:           int   *eptr;       // The offsets of each hyperedge
1763:           int   *eind;       // The vertices in each hyperedge, indexed by eptr
1764:           int   *hewgts;     // The hyperedge weights
1765:           int    nparts;     // The number of partitions
1766:           int    ubfactor;   // The allowed load imbalance (1-50)
1767:           int    options[9]; // Options
1768:           // Outputs
1769:           int    edgeCut;    // The number of edges cut by the partition
1770:           const Obj<ALE::Mesh::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);

1772:           if (topology->commRank() == 0) {
1773:             nvtxs      = bundle->heightStratum(1)->size();
1774:             vwgts      = NULL;
1775:             hewgts     = NULL;
1776:             nparts     = bundle->commSize();
1777:             ubfactor   = 5;
1778:             options[0] = 1;  // Use all defaults
1779:             options[1] = 10; // Number of bisections tested
1780:             options[2] = 1;  // Vertex grouping scheme
1781:             options[3] = 1;  // Objective function
1782:             options[4] = 1;  // V-cycle refinement
1783:             options[5] = 0;
1784:             options[6] = 0;
1785:             options[7] = 1; // Random seed
1786:             options[8] = 24; // Debugging level
1787:             assignment = new part_type[nvtxs];

1789:             if (bundle->commSize() == 1) {
1790:               PetscMemzero(assignment, nvtxs * sizeof(part_type));
1791:             } else {
1792:               ALE::New::Partitioner<bundle_type>::buildFaceCSR(bundle, dim, fNumbering, &nhedges, &eptr, &eind);
1793:               HMETIS_PartKway(nvtxs, nhedges, vwgts, eptr, eind, hewgts, nparts, ubfactor, options, assignment, &edgeCut);

1795:               delete [] eptr;
1796:               delete [] eind;
1797:             }
1798:             if (bundle->debug()) {for (int i = 0; i<nvtxs; i++) printf("[%d] %d\n", PetscGlobalRank, assignment[i]);}
1799:           } else {
1800:             assignment = NULL;
1801:           }
1802:           return assignment;
1803: #else
1804:           throw ALE::Exception("hmetis partitioner is not available.");
1805: #endif
1806:         };
1807:       };
1808:     };
1809: #endif
1810: #ifdef PETSC_HAVE_ZOLTAN
1811:     namespace Zoltan {
1812:       template<typename Bundle_>
1813:       class Partitioner {
1814:       public:
1815:         typedef Bundle_                          bundle_type;
1816:         typedef typename bundle_type::sieve_type sieve_type;
1817:         typedef typename bundle_type::point_type point_type;
1818:         typedef int                              part_type;
1819:       public:
1820:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1821:           throw ALE::Exception("Zoltan partition by cells not implemented");
1822:         };
1825:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1826:           // Outputs
1827:           float         version;           // The library version
1828:           int           changed;           // Did the partition change?
1829:           int           numGidEntries;     // Number of array entries for a single global ID (1)
1830:           int           numLidEntries;     // Number of array entries for a single local ID (1)
1831:           int           numImport;         // The number of imported points
1832:           ZOLTAN_ID_PTR import_global_ids; // The imported points
1833:           ZOLTAN_ID_PTR import_local_ids;  // The imported points
1834:           int          *import_procs;      // The proc each point was imported from
1835:           int          *import_to_part;    // The partition of each imported point
1836:           int           numExport;         // The number of exported points
1837:           ZOLTAN_ID_PTR export_global_ids; // The exported points
1838:           ZOLTAN_ID_PTR export_local_ids;  // The exported points
1839:           int          *export_procs;      // The proc each point was exported to
1840:           int          *export_to_part;    // The partition assignment of all local points
1841:           int          *assignment;        // The partition assignment of all local points
1842:           const Obj<typename bundle_type::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);

1844:           if (bundle->commSize() == 1) {
1845:             PetscMemzero(assignment, bundle->heightStratum(1)->size() * sizeof(part_type));
1846:           } else {
1847:             if (bundle->commRank() == 0) {
1848:               nvtxs_Zoltan = bundle->heightStratum(1)->size();
1849:               ALE::New::Partitioner<bundle_type>::buildFaceCSR(bundle, dim, fNumbering, &nhedges_Zoltan, &eptr_Zoltan, &eind_Zoltan);
1850:               assignment = new int[nvtxs_Zoltan];
1851:             } else {
1852:               nvtxs_Zoltan   = bundle->heightStratum(1)->size();
1853:               nhedges_Zoltan = 0;
1854:               eptr_Zoltan    = new int[1];
1855:               eind_Zoltan    = new int[1];
1856:               eptr_Zoltan[0] = 0;
1857:               assignment     = NULL;
1858:             }

1860:             int Zoltan_Initialize(0, NULL, &version);
1861:             struct Zoltan_Struct *zz = Zoltan_Create(bundle->comm());
1862:             // General parameters
1863:             Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
1864:             Zoltan_Set_Param(zz, "LB_METHOD", "PHG");
1865:             Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTITION");
1866:             // PHG parameters
1867:             Zoltan_Set_Param(zz, "PHG_OUTPUT_LEVEL", "2");
1868:             Zoltan_Set_Param(zz, "PHG_EDGE_SIZE_THRESHOLD", "1.0"); // Do not throw out dense edges
1869:             // Call backs
1870:             Zoltan_Set_Num_Obj_Fn(zz, getNumVertices_Zoltan, NULL);
1871:             Zoltan_Set_Obj_List_Fn(zz, getLocalElements_Zoltan, NULL);
1872:             Zoltan_Set_HG_Size_CS_Fn(zz, getHgSizes_Zoltan, NULL);
1873:             Zoltan_Set_HG_CS_Fn(zz, getHg_Zoltan, NULL);
1874:             // Debugging
1875:             //Zoltan_Generate_Files(zz, "zoltan.debug", 1, 0, 0, 1); // if using hypergraph callbacks

1877:             Zoltan_LB_Partition(zz, &changed, &numGidEntries, &numLidEntries,
1878:                                        &numImport, &import_global_ids, &import_local_ids, &import_procs, &import_to_part,
1879:                                        &numExport, &export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1880:             for(int v = 0; v < nvtxs_Zoltan; ++v) {
1881:               assignment[v] = export_to_part[v];
1882:             }
1883:             Zoltan_LB_Free_Part(&import_global_ids, &import_local_ids, &import_procs, &import_to_part);
1884:             Zoltan_LB_Free_Part(&export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1885:             Zoltan_Destroy(&zz);

1887:             delete [] eptr_Zoltan;
1888:             delete [] eind_Zoltan;
1889:           }
1890:           if (assignment) {for (int i=0; i<nvtxs_Zoltan; i++) printf("[%d] %d\n",PetscGlobalRank,assignment[i]);}
1891:           return assignment;
1892:         };
1893:       };
1894:     };
1895: #endif
1896:   }
1897: }
1898: #endif