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 <sieve/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:       std::map<point_type, point_type> numbers;
338:     protected:
339:       void createCells(const int height) {
340:         const Obj<typename Mesh::label_sequence>&     mcells   = mesh->heightStratum(height);
341:         const typename Mesh::label_sequence::iterator cEnd     = mcells->end();
342:         const int                                     numCells = mcells->size();
343:         int                                           c        = 0;

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

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

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

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

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

441:         if (size == 2) {
442:           offsets[cell+1]++;
443:         } else if ((size == 1) && (ovSize == 1)) {
444:           offsets[cell+1]++;
445:         }
446:       };
447:     };
448:     template<typename Sieve>
449:     class AdjVisitor {
450:     protected:
451:       typename Sieve::point_type cell;
452:       int                       *adjacency;
453:       const int                  cellOffset;
454:       int                        offset;
455:     public:
456:       AdjVisitor(int adj[], const bool zeroBase) : adjacency(adj), cellOffset(zeroBase ? 0 : 1), offset(0) {};
457:       void visitPoint(const typename Sieve::point_type& point) {};
458:       void visitArrow(const typename Sieve::arrow_type& arrow) {
459:         const int neighbor = arrow.target;

461:         if (neighbor != this->cell) {
462:           //std::cout << "Adding dual edge from " << cell << " to " << neighbor << std::endl;
463:           this->adjacency[this->offset++] = neighbor + this->cellOffset;
464:         }
465:       };
466:     public:
467:       void setCell(const typename Sieve::point_type cell) {this->cell = cell;};
468:       int  getOffset() {return this->offset;}
469:     };
470:     template<typename Mesh, int maxSize = 10>
471:     class FaceRecognizer {
472:     public:
473:       typedef typename Mesh::point_type point_type;
474:     protected:
475:       int numCases;
476:       int numFaceVertices[maxSize];
477:     public:
478:       FaceRecognizer(Mesh& mesh, const int debug = 0) : numCases(0) {
479:         if (mesh.depth() != 1) {throw ALE::Exception("Only works for depth 1 meshes");}
480:         const int                                 dim   = mesh.getDimension();
481:         const Obj<typename Mesh::sieve_type>&     sieve = mesh.getSieve();
482:         const Obj<typename Mesh::label_sequence>& cells = mesh.heightStratum(0);
483:         std::set<int>                             cornersSeen;

485:         if (debug) {std::cout << "Building Recognizer" << std::endl;}
486:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
487:           const int corners = sieve->getConeSize(*c_iter);

489:           if (cornersSeen.find(corners) == cornersSeen.end()) {
490:             if (numCases >= maxSize) {throw ALE::Exception("Exceeded maximum number of cases");}
491:             cornersSeen.insert(corners);

493:             if (corners == dim+1) {
494:               numFaceVertices[numCases] = dim;
495:               if (debug) {std::cout << "  Recognizing simplices" << std::endl;}
496:             } else if ((dim == 1) && (corners == 3)) {
497:               numFaceVertices[numCases] = 3;
498:               if (debug) {std::cout << "  Recognizing quadratic edges" << std::endl;}
499:             } else if ((dim == 2) && (corners == 4)) {
500:               numFaceVertices[numCases] = 2;
501:               if (debug) {std::cout << "  Recognizing quads" << std::endl;}
502:             } else if ((dim == 2) && (corners == 6)) {
503:               numFaceVertices[numCases] = 3;
504:               if (debug) {std::cout << "  Recognizing tri and quad cohesive Lagrange cells" << std::endl;}
505:             } else if ((dim == 2) && (corners == 9)) {
506:               numFaceVertices[numCases] = 3;
507:               if (debug) {std::cout << "  Recognizing quadratic quads and quadratic quad cohesive Lagrange cells" << std::endl;}
508:             } else if ((dim == 3) && (corners == 6)) {
509:               numFaceVertices[numCases] = 4;
510:               if (debug) {std::cout << "  Recognizing tet cohesive cells" << std::endl;}
511:             } else if ((dim == 3) && (corners == 8)) {
512:               numFaceVertices[numCases] = 4;
513:               if (debug) {std::cout << "  Recognizing hexes" << std::endl;}
514:             } else if ((dim == 3) && (corners == 9)) {
515:               numFaceVertices[numCases] = 6;
516:               if (debug) {std::cout << "  Recognizing tet cohesive Lagrange cells" << std::endl;}
517:             } else if ((dim == 3) && (corners == 10)) {
518:               numFaceVertices[numCases] = 6;
519:               if (debug) {std::cout << "  Recognizing quadratic tets" << std::endl;}
520:             } else if ((dim == 3) && (corners == 12)) {
521:               numFaceVertices[numCases] = 6;
522:               if (debug) {std::cout << "  Recognizing hex cohesive Lagrange cells" << std::endl;}
523:             } else if ((dim == 3) && (corners == 18)) {
524:               numFaceVertices[numCases] = 6;
525:               if (debug) {std::cout << "  Recognizing quadratic tet cohesive Lagrange cells" << std::endl;}
526:             } else if ((dim == 3) && (corners == 27)) {
527:               numFaceVertices[numCases] = 9;
528:               if (debug) {std::cout << "  Recognizing quadratic hexes and quadratic hex cohesive Lagrange cells" << std::endl;}
529:             } else {
530:               throw ALE::Exception("Could not determine number of face vertices");
531:             }
532:             ++numCases;
533:           }
534:         }
535:       };
536:       ~FaceRecognizer() {};
537:     public:
538:       bool operator()(point_type cellA, point_type cellB, int meetSize) {
539:         // Could concievably make this depend on the cells, but it seems slow
540:         for(int i = 0; i < numCases; ++i) {
541:           if (meetSize == numFaceVertices[i]) {
542:             //std::cout << "  Recognized case " << i <<"("<<numFaceVertices[i]<<") for <" << cellA <<","<< cellB << ">" << std::endl;
543:             return true;
544:           }
545:         }
546:         return false;
547:       };
548:     };
549:     template<typename Sieve, typename Manager, typename Recognizer>
550:     class MeetVisitor {
551:     public:
552:       typedef std::set<typename Sieve::point_type> neighbors_type;
553:     protected:
554:       const Sieve& sieve;
555:       Manager& manager;
556:       Recognizer& faceRecognizer;
557:       const int numCells;
558:       neighbors_type *neighborCells;
559:       typename ISieveVisitor::PointRetriever<Sieve> *pR;
560:       typename Sieve::point_type cell;
561:     public:
562:       MeetVisitor(const Sieve& s, Manager& manager, Recognizer& faceRecognizer, const int n) : sieve(s), manager(manager), faceRecognizer(faceRecognizer), numCells(n) {
563:         this->neighborCells = new std::set<typename Sieve::point_type>[numCells];
564:         this->pR            = new typename ISieveVisitor::PointRetriever<Sieve>(this->sieve.getMaxConeSize());
565:       };
566:       ~MeetVisitor() {delete [] this->neighborCells; delete this->pR;};
567:       void visitArrow(const typename Sieve::arrow_type& arrow) {};
568:       void visitPoint(const typename Sieve::point_type& point) {
569:         const typename Sieve::point_type& neighbor = point;

571:         if (this->cell == neighbor) return;
572:         this->pR->clear();
573:         this->sieve.meet(this->cell, neighbor, *this->pR);
574:         if (faceRecognizer(this->cell, neighbor, this->pR->getSize())) {
575:           this->neighborCells[this->manager.getNumber(this->cell)].insert(this->manager.getNumber(neighbor));
576:         }
577:       };
578:     public:
579:       void setCell(const typename Sieve::point_type& c) {this->cell = c;};
580:       const neighbors_type *getNeighbors() {return this->neighborCells;};
581:     };
582:   public: // Creating overlaps
583:     // Create a partition point overlap for distribution
584:     //   This is the default overlap which comes from distributing a serial mesh on process 0
585:     template<typename SendOverlap, typename RecvOverlap>
586:     static void createDistributionPartOverlap(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
587:       const int rank = sendOverlap->commRank();
588:       const int size = sendOverlap->commSize();

590:       if (rank == 0) {
591:         for(int p = 1; p < size; p++) {
592:           // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
593:           sendOverlap->addCone(p, p, p);
594:         }
595:       }
596:       if (rank != 0) {
597:         // The arrow is from remote partition point rank (color) on rank 0 (source) to local partition point rank (target)
598:         recvOverlap->addCone(0, rank, rank);
599:       }
600:     }
601:     // Create a mesh point overlap for distribution
602:     //   A local numbering is created for the remote points
603:     //   This is the default overlap which comes from distributing a serial mesh on process 0
604:     template<typename Section, typename RecvPartOverlap, typename Renumbering, typename SendOverlap, typename RecvOverlap>
605:     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) {
606:       const typename Section::chart_type& chart = partition->getChart();
607:       int numRanks = 0;

609:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
610:         if (*p_iter != sendOverlap->commRank() && partition->getFiberDimension(*p_iter)) numRanks++;
611:       }
612:       sendOverlap->setBaseSize(numRanks); // setNumRanks
613:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
614:         const int coneSize = partition->getFiberDimension(*p_iter);

616:         if (*p_iter != sendOverlap->commRank() && coneSize) {
617:           sendOverlap->setConeSize(*p_iter, coneSize); // setNumPoints
618:         }
619:       }
620:       sendOverlap->assemble();
621:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
622:         if (*p_iter == sendOverlap->commRank()) continue;
623:         const typename Section::value_type *points     = partition->restrictPoint(*p_iter);
624:         const int                           numPoints  = partition->getFiberDimension(*p_iter);

626:         for(int i = 0; i < numPoints; ++i) {
627:           // Notice here that we do not know the local renumbering (but we do not use it)
628:           sendOverlap->addArrow(points[i], *p_iter, points[i]);
629:         }
630:       }
631:       if (sendOverlap->debug()) {sendOverlap->view("Send mesh overlap");}
632:       const typename RecvPartOverlap::capSequence::iterator rBegin = recvPartOverlap->capBegin();
633:       const typename RecvPartOverlap::capSequence::iterator rEnd   = recvPartOverlap->capEnd();

635:       recvOverlap->setCapSize(recvPartOverlap->getCapSize()); // setNumRanks
636:       for(typename RecvPartOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
637:         const int                                                 rank   = *r_iter;
638:         const typename RecvPartOverlap::supportSequence::iterator pBegin = recvPartOverlap->supportBegin(*r_iter);
639:         const typename RecvPartOverlap::supportSequence::iterator pEnd   = recvPartOverlap->supportEnd(*r_iter);

641:         for(typename RecvPartOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
642:           const typename Section::point_type& remotePartPoint = p_iter.color();
643:           const int                           numPoints       = overlapPartition->getFiberDimension(remotePartPoint);

645:           recvOverlap->setSupportSize(rank, numPoints); // setNumPoints
646:         }
647:       }
648:       recvOverlap->assemble();

650:       for(typename RecvPartOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
651:         const int                                                 rank   = *r_iter;
652:         const typename RecvPartOverlap::supportSequence::iterator pBegin = recvPartOverlap->supportBegin(*r_iter);
653:         const typename RecvPartOverlap::supportSequence::iterator pEnd   = recvPartOverlap->supportEnd(*r_iter);

655:         for(typename RecvPartOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
656:           const typename Section::point_type& remotePartPoint = p_iter.color();
657:           const typename Section::value_type *points          = overlapPartition->restrictPoint(remotePartPoint);
658:           const int                           numPoints       = overlapPartition->getFiberDimension(remotePartPoint);

660:           for(int i = 0; i < numPoints; ++i) {
661:             recvOverlap->addArrow(rank, renumbering[points[i]], points[i]);
662:           }
663:         }
664:       }
665:       if (recvOverlap->debug()) {recvOverlap->view("Receive mesh overlap");}
666:     }
667:     // Create a partition point overlap from a partition
668:     //   The intention is to create an overlap which enables exchange of redistribution information
669:     template<typename Section, typename SendOverlap, typename RecvOverlap>
670:     static void createPartitionPartOverlap(const Obj<Section>& partition, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
671:       const typename Section::chart_type& chart      = partition->getChart();
672:       const int                           rank       = partition->commRank();
673:       const int                           size       = partition->commSize();
674:       int                                *adj        = new int[size];
675:       int                                *recvCounts = new int[size];
676:       int                                 numNeighbors;

678:       for(int p = 0; p < size; ++p) {
679:         adj[p]        = 0;
680:         recvCounts[p] = 1;
681:       }
682:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart->end(); ++p_iter) {
683:         const typename Section::value_type& p = partition->restrictPoint(*p_iter)[0];
684:         // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
685:         sendOverlap->addCone(p, p, p);
686:         adj[p] = 1;
687:       }
688:       MPI_Reduce_Scatter(adj, &numNeighbors, recvCounts, size, MPI_INT, MPI_SUM, partition->comm());
689:       MPI_Request *recvRequests = new MPI_Request[numNeighbors];
690:       int          dummy        = 0;

692:       // TODO: Get a unique tag
693:       for(int n = 0; n < numNeighbors; ++n) {
694:         MPI_Irecv(&dummy, 1, MPI_INT, MPI_ANY_SOURCE, 1, partition->comm(), &recvRequests[n]);
695:       }
696:       const Obj<typename SendOverlap::traits::baseSequence>      ranks        = sendOverlap->base();
697:       const typename SendOverlap::traits::baseSequence::iterator rEnd         = ranks->end();
698:       MPI_Request                                               *sendRequests = new MPI_Request[ranks->size()];
699:       int                                                        s            = 0;

701:       for(typename SendOverlap::traits::baseSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter, ++s) {
702:         MPI_Isend(&dummy, 1, MPI_INT, *r_iter, 1, partition->comm(), &sendRequests[s]);
703:       }
704:       for(int n = 0; n < numNeighbors; ++n) {
705:         MPI_Status status;
706:         int        idx;

708:         MPI_Waitany(numNeighbors, recvRequests, &idx, &status);
709:         // The arrow is from remote partition point rank (color) on rank p (source) to local partition point rank (target)
710:         recvOverlap->addCone(status.MPI_SOURCE, rank, rank);
711:       }
712:       MPI_Waitall(ranks->size(), sendRequests, MPI_STATUSES_IGNORE);
713:       delete [] sendRequests;
714:       delete [] recvRequests;
715:       delete [] adj;
716:       delete [] recvCounts;
717:     }
718:   public: // Building CSR meshes
719:     // This produces the dual graph (each cell is a vertex and each face is an edge)
720:     //   numbering:   A contiguous numbering of the cells (not yet included)
721:     //   numVertices: The number of vertices in the graph (cells in the mesh)
722:     //   adjacency:   The vertices adjacent to each vertex (cells adjacent to each mesh cell)
723:     // - We allow an exception to contiguous numbering.
724:     //   If the cell id > numElements, we assign a new number starting at
725:     //     the top and going downward. I know these might not match up with
726:     //     the iterator order, but we can fix it later.
727:     template<typename Mesh>
728:     static void buildDualCSR(const Obj<Mesh>& mesh, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
729:       const Obj<typename Mesh::sieve_type>&         sieve        = mesh->getSieve();
730:       const Obj<typename Mesh::label_sequence>&     cells        = mesh->heightStratum(0);
731:       const typename Mesh::label_sequence::iterator cEnd         = cells->end();
732:       const int                                     numCells     = cells->size();
733:       int                                           newCell      = numCells;
734:       Obj<typename Mesh::sieve_type>                overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
735:       int                                           offset       = 0;
736:       const int                                     cellOffset   = zeroBase ? 0 : 1;
737:       const int                                     dim          = mesh->getDimension();
738:       std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;

740:       // TODO: This is necessary for parallel partitioning
741:       //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
742:       if (numCells == 0) {
743:         *numVertices = 0;
744:         *offsets     = NULL;
745:         *adjacency   = NULL;
746:         return;
747:       }
748:       int *off = alloc_type().allocate(numCells+1);
749:       int *adj;
750:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
751:       if (mesh->depth() == dim) {
752:         int c = 1;

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

758:           off[c] = off[c-1];
759:           for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
760:             if (sieve->support(*f_iter)->size() == 2) {
761:               off[c]++;
762:             } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
763:               off[c]++;
764:             }
765:           }
766:         }
767:         adj = alloc_type().allocate(off[numCells]);
768:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
769:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
770:           const Obj<typename Mesh::sieve_type::traits::coneSequence>&     faces = sieve->cone(*c_iter);
771:           const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd  = faces->end();

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

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

783:             for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = oNeighbors->begin(); n_iter != onEnd; ++n_iter) {
784:               adj[offset++] = *n_iter + cellOffset;
785:             }
786:           }
787:         }
788:       } else if (mesh->depth() == 1) {
789:         std::set<typename Mesh::point_type> *neighborCells = new std::set<typename Mesh::point_type>[numCells];
790:         const int                            corners       = sieve->cone(*cells->begin())->size();
791:         int                                  faceVertices;

793:         if (corners == dim+1) {
794:           faceVertices = dim;
795:         } else if ((dim == 2) && (corners == 4)) {
796:           faceVertices = 2;
797:         } else if ((dim == 3) && (corners == 8)) {
798:           faceVertices = 4;
799:         } else if ((dim == 2) && (corners == 6)) {
800:           faceVertices = 3;
801:         } else if ((dim == 2) && (corners == 9)) {
802:           faceVertices = 3;
803:         } else if ((dim == 3) && (corners == 10)) {
804:           faceVertices = 6;
805:         } else if ((dim == 3) && (corners == 27)) {
806:           faceVertices = 9;
807:         } else {
808:           throw ALE::Exception("Could not determine number of face vertices");
809:         }
810:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
811:             const Obj<typename Mesh::sieve_type::traits::coneSequence>&     vertices = sieve->cone(*c_iter);
812:             const typename Mesh::sieve_type::traits::coneSequence::iterator vEnd     = vertices->end();

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

818:               for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
819:                 if (*c_iter == *n_iter) continue;
820:                 if ((int) sieve->nMeet(*c_iter, *n_iter, 1)->size() == faceVertices) {
821:                   if ((*c_iter < numCells) && (*n_iter < numCells)) {
822:                     neighborCells[*c_iter].insert(*n_iter);
823:                   } else {
824:                     typename Mesh::point_type e = *c_iter, n = *n_iter;

826:                     if (*c_iter >= numCells) {
827:                       if (newCells.find(*c_iter) == newCells.end()) newCells[*c_iter] = --newCell;
828:                       e = newCells[*c_iter];
829:                     }
830:                     if (*n_iter >= numCells) {
831:                       if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
832:                       n = newCells[*n_iter];
833:                     }
834:                     neighborCells[e].insert(n);
835:                   }
836:                 }
837:               }
838:             }
839:           }
840:           off[0] = 0;
841:           for(int c = 1; c <= numCells; c++) {
842:             off[c] = neighborCells[c-1].size() + off[c-1];
843:           }
844:           adj = alloc_type().allocate(off[numCells]);
845:           for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
846:           for(int c = 0; c < numCells; c++) {
847:             for(typename std::set<typename Mesh::point_type>::iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
848:               adj[offset++] = *n_iter + cellOffset;
849:             }
850:           }
851:           delete [] neighborCells;
852:       } else {
853:         throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
854:       }
855:       if (offset != off[numCells]) {
856:         ostringstream msg;
857:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
858:         throw ALE::Exception(msg.str().c_str());
859:       }
860:       *numVertices = numCells;
861:       *offsets     = off;
862:       *adjacency   = adj;
863:     }
864:     template<typename Mesh>
865:     static void buildDualCSRV(const Obj<Mesh>& mesh, MeshManager<Mesh>& manager, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
866:       const Obj<typename Mesh::sieve_type>&         sieve        = mesh->getSieve();
867:       const Obj<typename Mesh::label_sequence>&     cells        = mesh->heightStratum(0);
868:       const typename Mesh::label_sequence::iterator cEnd         = cells->end();
869:       const int                                     numCells     = cells->size();
870:       Obj<typename Mesh::sieve_type>                overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
871:       int                                           offset       = 0;
872:       const int                                     cellOffset   = zeroBase ? 0 : 1;
873:       const int                                     dim          = mesh->getDimension();
874:       std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;

876:       // TODO: This is necessary for parallel partitioning
877:       //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
878:       overlapSieve->setChart(sieve->getChart());
879:       overlapSieve->allocate();
880:       if (numCells == 0) {
881:         *numVertices = 0;
882:         *offsets     = NULL;
883:         *adjacency   = NULL;
884:         return;
885:       }
886:       int *off = alloc_type().allocate(numCells+1);
887:       int *adj;
888:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
889:       if (mesh->depth() == dim) {
890:         OffsetVisitor<typename Mesh::sieve_type> oV(*sieve, *overlapSieve, off);

892:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
893:           sieve->cone(*c_iter, oV);
894:         }
895:         for(int p = 1; p <= numCells; ++p) {
896:           off[p] = off[p] + off[p-1];
897:         }
898:         adj = alloc_type().allocate(off[numCells]);
899:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
900:         AdjVisitor<typename Mesh::sieve_type> aV(adj, zeroBase);
901:         ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > sV(*sieve, aV);
902:         ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > ovSV(*overlapSieve, aV);

904:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
905:           aV.setCell(*c_iter);
906:           sieve->cone(*c_iter, sV);
907:           sieve->cone(*c_iter, ovSV);
908:         }
909:         offset = aV.getOffset();
910:       } else if (mesh->depth() == 1) {
911:         typedef MeetVisitor<typename Mesh::sieve_type, MeshManager<Mesh>, FaceRecognizer<Mesh> > mv_type;
912:         typedef typename ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, mv_type>       sv_type;
913:         FaceRecognizer<Mesh> faceRecognizer(*mesh);

915:         mv_type mV(*sieve, manager, faceRecognizer, numCells);
916:         sv_type sV(*sieve, mV);

918:         for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
919:           mV.setCell(*c_iter);
920:           sieve->cone(*c_iter, sV);
921:         }
922:         const typename mv_type::neighbors_type *neighborCells = mV.getNeighbors();

924:         off[0] = 0;
925:         for(int c = 1; c <= numCells; c++) {
926:           off[c] = neighborCells[c-1].size() + off[c-1];
927:         }
928:         adj = alloc_type().allocate(off[numCells]);
929:         for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
930:         for(int c = 0; c < numCells; c++) {
931:           for(typename mv_type::neighbors_type::const_iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
932:             //std::cout << "Adding dual edge from " << c << " to " << *n_iter << std::endl;
933:             adj[offset++] = *n_iter + cellOffset;
934:           }
935:         }
936:       } else {
937:         throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
938:       }
939:       if (offset != off[numCells]) {
940:         ostringstream msg;
941:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
942:         throw ALE::Exception(msg.str().c_str());
943:       }
944:       *numVertices = numCells;
945:       *offsets     = off;
946:       *adjacency   = adj;
947:     }
948:     // This produces a hypergraph (each face is a vertex and each cell is a hyperedge)
949:     //   numbering: A contiguous numbering of the faces
950:     //   numEdges:  The number of edges in the hypergraph
951:     //   adjacency: The vertices in each edge
952:     template<typename Mesh>
953:     static void buildFaceDualCSR(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
954:       const Obj<typename Mesh::sieve_type>&         sieve = mesh->getSieve();
955:       const Obj<typename Mesh::label_sequence>&     cells = mesh->heightStratum(0);
956:       const typename Mesh::label_sequence::iterator cEnd  = cells->end();
957:       const int faceOffset = zeroBase ? 0 : 1;
958:       int       numCells = cells->size();
959:       int       c        = 1;

961:       if (mesh->depth() != mesh->getDimension()) {throw ALE::Exception("Not yet implemented for non-interpolated meshes");}
962:       int *off = alloc_type().allocate(numCells+1);
963:       for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
964:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
965:         off[c] = sieve->cone(*c_iter)->size() + off[c-1];
966:       }
967:       int *adj = alloc_type().allocate(off[numCells]);
968:       for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
969:       int  offset = 0;
970:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
971:         const Obj<typename Mesh::sieve_type::traits::coneSequence>&     faces = sieve->cone(*c_iter);
972:         const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd  = faces->end();

974:         for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter, ++offset) {
975:           adj[offset] = numbering->getIndex(*f_iter) + faceOffset;
976:         }
977:       }
978:       if (offset != off[numCells]) {
979:         ostringstream msg;
980:         msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
981:         throw ALE::Exception(msg.str().c_str());
982:       }
983:       *numEdges  = numCells;
984:       *offsets   = off;
985:       *adjacency = adj;
986:     }
987:     template<typename Mesh>
988:     static void buildFaceDualCSRV(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
989:       throw ALE::Exception("Not implemented");
990:     }
991:     static void destroyCSR(int numPoints, int *offsets, int *adjacency) {
992:       if (adjacency) {
993:         for(int i = 0; i < offsets[numPoints]; ++i) {alloc_type().destroy(adjacency+i);}
994:         alloc_type().deallocate(adjacency, offsets[numPoints]);
995:       }
996:       if (offsets) {
997:         for(int i = 0; i < numPoints+1; ++i) {alloc_type().destroy(offsets+i);}
998:         alloc_type().deallocate(offsets, numPoints+1);
999:       }
1000:     }
1001:     template<typename OldSection, typename Partition, typename Renumbering, typename NewSection>
1002:     static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<NewSection>& newSection) {
1003:       const typename Partition::value_type *points    = partition->restrictPoint(oldSection->commRank());
1004:       const int                             numPoints = partition->getFiberDimension(oldSection->commRank());

1006:       for(int p = 0; p < numPoints; ++p) {
1007:         if (oldSection->hasPoint(points[p])) {
1008:           newSection->setFiberDimension(renumbering[points[p]], oldSection->getFiberDimension(points[p]));
1009:         }
1010:       }
1011:       if (numPoints) {newSection->allocatePoint();}
1012:       for(int p = 0; p < numPoints; ++p) {
1013:         if (oldSection->hasPoint(points[p])) {
1014:           newSection->updatePointAll(renumbering[points[p]], oldSection->restrictPoint(points[p]));
1015:         }
1016:       }
1017:     }
1018:     // Specialize to ArrowSection
1019:     template<typename OldSection, typename Partition, typename Renumbering>
1020:     static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<UniformSection<MinimalArrow<int,int>,int> >& newSection) {
1021:       typedef UniformSection<MinimalArrow<int,int>,int> NewSection;
1022:       const typename Partition::value_type    *points    = partition->restrictPoint(oldSection->commRank());
1023:       const int                                numPoints = partition->getFiberDimension(oldSection->commRank());
1024:       const typename OldSection::chart_type&   oldChart  = oldSection->getChart();
1025:       std::set<typename Partition::value_type> myPoints;

1027:       for(int p = 0; p < numPoints; ++p) {
1028:         myPoints.insert(points[p]);
1029:       }
1030:       for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
1031:         if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
1032:           newSection->setFiberDimension(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), oldSection->getFiberDimension(*c_iter));
1033:         }
1034:       }
1035:       if (oldChart.size()) {newSection->allocatePoint();}
1036:       for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
1037:         if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
1038:           const typename OldSection::value_type *values = oldSection->restrictPoint(*c_iter);

1040:           newSection->updatePointAll(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), values);
1041:         }
1042:       }
1043:     }
1044:     template<typename Sifter, typename Section, typename Renumbering>
1045:     static void createLocalSifter(const Obj<Sifter>& sifter, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sifter>& localSifter) {
1046:       const typename Section::value_type *points    = partition->restrictPoint(sifter->commRank());
1047:       const int                           numPoints = partition->getFiberDimension(sifter->commRank());

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

1053:         for(typename Sifter::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
1054:           localSifter->addArrow(*c_iter, renumbering[points[p]]);
1055:         }
1056:       }
1057:     }
1058:     template<typename Sieve, typename Section, typename Renumbering>
1059:     static void createLocalSieve(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1060:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1061:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());

1063:       for(int p = 0; p < numPoints; ++p) {
1064:         Obj<typename Sieve::coneSet> current = new typename Sieve::coneSet();
1065:         Obj<typename Sieve::coneSet> next    = new typename Sieve::coneSet();
1066:         Obj<typename Sieve::coneSet> tmp;

1068:         current->insert(points[p]);
1069:         while(current->size()) {
1070:           for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1071:             const Obj<typename Sieve::traits::coneSequence>& cone = sieve->cone(*p_iter);
1072: 
1073:             for(typename Sieve::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1074:               localSieve->addArrow(renumbering[*c_iter], renumbering[*p_iter], c_iter.color());
1075:               next->insert(*c_iter);
1076:             }
1077:           }
1078:           tmp = current; current = next; next = tmp;
1079:           next->clear();
1080:         }
1081:         if (height) {
1082:           current->insert(points[p]);
1083:           while(current->size()) {
1084:             for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1085:               const Obj<typename Sieve::traits::supportSequence>& support = sieve->support(*p_iter);
1086: 
1087:               for(typename Sieve::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1088:                 localSieve->addArrow(renumbering[*p_iter], renumbering[*s_iter], s_iter.color());
1089:                 next->insert(*s_iter);
1090:               }
1091:             }
1092:             tmp = current; current = next; next = tmp;
1093:             next->clear();
1094:           }
1095:         }
1096:       }
1097:     }
1098:     template<typename Mesh, typename Section, typename Renumbering>
1099:     static void createLocalMesh(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1100:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
1101:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1103:       createLocalSieve(sieve, partition, renumbering, localSieve, height);
1104:     }
1105:     template<typename Sieve, typename Section, typename Renumbering>
1106:     static void sizeLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1107:       typedef std::set<typename Sieve::point_type> pointSet;
1108:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1109:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1110:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1111:       const pointSet                      pSet(points, &points[numPoints]);
1112:       ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> fV(pSet, renumbering, maxSize);

1114:       for(int p = 0; p < numPoints; ++p) {
1115:         sieve->cone(points[p], fV);
1116:         localSieve->setConeSize(renumbering[points[p]], fV.getSize());
1117:         fV.clear();
1118:         sieve->support(points[p], fV);
1119:         localSieve->setSupportSize(renumbering[points[p]], fV.getSize());
1120:         fV.clear();
1121:       }
1122:     }
1123:     template<typename Mesh, typename Section, typename Renumbering>
1124:     static void sizeLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1125:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
1126:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1128:       sizeLocalSieveV(sieve, partition, renumbering, localSieve, height);
1129:     }
1130:     template<typename Sieve, typename Section, typename Renumbering>
1131:     static void createLocalLabelV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1132:       typedef std::set<typename Sieve::point_type> pointSet;
1133:       typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1134:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1135:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1136:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1137:       typename Sieve::point_type         *oPoints   = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1138:       int                                *oOrients  = new int[std::max(1, sieve->getMaxConeSize())];
1139:       const pointSet                      pSet(points, &points[numPoints]);
1140:       visitor_type fV(pSet, renumbering, maxSize);

1142:       for(int p = 0; p < numPoints; ++p) {
1143:         fV.useRenumbering(false);
1144:         sieve->orientedCone(points[p], fV);
1145:         const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1146:         const int                                         n = fV.getOrientedSize();
1147:         for(int i = 0; i < n; ++i) {
1148:           oPoints[i]  = q[i].first;
1149:           oOrients[i] = q[i].second;
1150:         }
1151:         localSieve->setCone(oPoints, renumbering[points[p]]);
1152:         localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1153:         fV.clear();
1154:         fV.useRenumbering(true);
1155:         sieve->support(points[p], fV);
1156:         if (fV.getSize()) {localSieve->setSupport(points[p], fV.getPoints());}
1157:         fV.clear();
1158:       }
1159:       delete [] oPoints;
1160:       delete [] oOrients;
1161:     }
1162:     template<typename Sieve, typename Section, typename Renumbering>
1163:     static void createLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1164:       typedef std::set<typename Sieve::point_type> pointSet;
1165:       typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1166:       const typename Section::value_type *points    = partition->restrictPoint(sieve->commRank());
1167:       const int                           numPoints = partition->getFiberDimension(sieve->commRank());
1168:       int                                 maxSize   = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1169:       typename Sieve::point_type         *oPoints   = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1170:       int                                *oOrients  = new int[std::max(1, sieve->getMaxConeSize())];
1171:       const pointSet                      pSet(points, &points[numPoints]);
1172:       visitor_type fV(pSet, renumbering, maxSize);

1174:       for(int p = 0; p < numPoints; ++p) {
1175:         ///sieve->cone(points[p], fV);
1176:         ///localSiaeve->setCone(fV.getPoints(), renumbering[points[p]]);
1177:         sieve->orientedCone(points[p], fV);
1178:         const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1179:         const int                                         n = fV.getOrientedSize();
1180:         for(int i = 0; i < n; ++i) {
1181:           oPoints[i]  = q[i].first;
1182:           oOrients[i] = q[i].second;
1183:         }
1184:         localSieve->setCone(oPoints, renumbering[points[p]]);
1185:         localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1186:         fV.clear();
1187:         sieve->support(points[p], fV);
1188:         localSieve->setSupport(renumbering[points[p]], fV.getPoints());
1189:         fV.clear();
1190:       }
1191:       delete [] oPoints;
1192:       delete [] oOrients;
1193:     }
1194:     template<typename Mesh, typename Section, typename Renumbering>
1195:     static void createLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1196:       const Obj<typename Mesh::sieve_type>& sieve      = mesh->getSieve();
1197:       const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();

1199:       createLocalSieveV(sieve, partition, renumbering, localSieve, height);
1200:     }
1201:   public: // Partitioning
1202:     //   partition:    Should be properly allocated on input
1203:     //   height:       Height of the point set to uniquely partition
1204:     // TODO: Could rebind assignment section to the type of the output
1205:     template<typename Mesh, typename Section>
1206:     static void createPartition(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1207:       MeshManager<Mesh> manager(mesh, height);
1208:       int              *start     = NULL;
1209:       int              *adjacency = NULL;

1211:       if (height == 0) {
1212:         int numVertices;

1214:         buildDualCSR(mesh, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1215:         GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1216:         destroyCSR(numVertices, start, adjacency);
1217:       } else if (height == 1) {
1218:         int numEdges;

1220:         buildFaceDualCSR(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1221:         GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1222:         destroyCSR(numEdges, start, adjacency);
1223:       } else {
1224:         throw ALE::Exception("Invalid partition height");
1225:       }
1226:     }
1227:     template<typename Mesh, typename Section>
1228:     static void createPartitionV(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1229:       MeshManager<Mesh> manager(mesh, height);
1230:       int              *start     = NULL;
1231:       int              *adjacency = NULL;

1233:       PETSc::Log::Event("PartitionCreate").begin();
1234:       if (height == 0) {
1235:         int numVertices;

1237:         buildDualCSRV(mesh, manager, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1238:         GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1239:         destroyCSR(numVertices, start, adjacency);
1240:       } else if (height == 1) {
1241:         int numEdges;

1243:         throw ALE::Exception("Not yet implemented");
1244: #if 0
1245:         buildFaceDualCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1246: #endif
1247:         GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1248:         destroyCSR(numEdges, start, adjacency);
1249:       } else {
1250:         throw ALE::Exception("Invalid partition height");
1251:       }
1252:       PETSc::Log::Event("PartitionCreate").end();
1253:     }
1254:     // Add in the points in the closure (and star) of the partitioned points
1255:     template<typename Mesh, typename Section>
1256:     static void createPartitionClosure(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1257:       const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1258:       const typename Section::chart_type&   chart = pointPartition->getChart();
1259:       size_t                                size  = 0;

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

1266:         // TODO: Use Quiver's closure() here instead
1267:         for(int p = 0; p < numPoints; ++p) {
1268:           Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1269:           Obj<typename Mesh::sieve_type::coneSet> next    = new typename Mesh::sieve_type::coneSet();
1270:           Obj<typename Mesh::sieve_type::coneSet> tmp;

1272:           current->insert(points[p]);
1273:           closure.insert(points[p]);
1274:           while(current->size()) {
1275:             for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1276:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1277: 
1278:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1279:                 closure.insert(*c_iter);
1280:                 next->insert(*c_iter);
1281:               }
1282:             }
1283:             tmp = current; current = next; next = tmp;
1284:             next->clear();
1285:           }
1286:           if (height) {
1287:             current->insert(points[p]);
1288:             while(current->size()) {
1289:               for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1290:                 const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1291: 
1292:                 for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1293:                   closure.insert(*s_iter);
1294:                   next->insert(*s_iter);
1295:                 }
1296:               }
1297:               tmp = current; current = next; next = tmp;
1298:               next->clear();
1299:             }
1300:           }
1301:         }
1302:         partition->setFiberDimension(*r_iter, closure.size());
1303:         size = std::max(size, closure.size());
1304:       }
1305:       partition->allocatePoint();
1306:       typename Section::value_type *values = new typename Section::value_type[size];

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

1313:         // TODO: Use Quiver's closure() here instead
1314:         for(int p = 0; p < numPoints; ++p) {
1315:           Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1316:           Obj<typename Mesh::sieve_type::coneSet> next    = new typename Mesh::sieve_type::coneSet();
1317:           Obj<typename Mesh::sieve_type::coneSet> tmp;

1319:           current->insert(points[p]);
1320:           closure.insert(points[p]);
1321:           while(current->size()) {
1322:             for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1323:               const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1324: 
1325:               for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1326:                 closure.insert(*c_iter);
1327:                 next->insert(*c_iter);
1328:               }
1329:             }
1330:             tmp = current; current = next; next = tmp;
1331:             next->clear();
1332:           }
1333:           if (height) {
1334:             current->insert(points[p]);
1335:             while(current->size()) {
1336:               for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1337:                 const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1338: 
1339:                 for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1340:                   closure.insert(*s_iter);
1341:                   next->insert(*s_iter);
1342:                 }
1343:               }
1344:               tmp = current; current = next; next = tmp;
1345:               next->clear();
1346:             }
1347:           }
1348:         }
1349:         int i = 0;

1351:         for(typename std::set<typename Section::value_type>::const_iterator p_iter = closure.begin(); p_iter != closure.end(); ++p_iter, ++i) {
1352:           values[i] = *p_iter;
1353:         }
1354:         partition->updatePoint(*r_iter, values);
1355:       }
1356:       delete [] values;
1357:     }
1358:     template<typename Mesh, typename Section>
1359:     static void createPartitionClosureV(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1360:       typedef ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type> visitor_type;
1361:       const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1362:       const typename Section::chart_type&   chart = pointPartition->getChart();
1363:       size_t                                size  = 0;

1365:       PETSc::Log::Event("PartitionClosure").begin();
1366:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1367:         const typename Section::value_type *points    = pointPartition->restrictPoint(*r_iter);
1368:         const int                           numPoints = pointPartition->getFiberDimension(*r_iter);
1369:         typename visitor_type::visitor_type nV;
1370:         visitor_type                        cV(*sieve, nV);

1372:         for(int p = 0; p < numPoints; ++p) {
1373:           sieve->cone(points[p], cV);
1374:           if (height) {
1375:             cV.setIsCone(false);
1376:             sieve->support(points[p], cV);
1377:           }
1378:         }
1379:         partition->setFiberDimension(*r_iter, cV.getPoints().size());
1380:         size = std::max(size, cV.getPoints().size());
1381:       }
1382:       partition->allocatePoint();
1383:       typename Section::value_type *values = new typename Section::value_type[size];

1385:       for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1386:         const typename Section::value_type *points    = pointPartition->restrictPoint(*r_iter);
1387:         const int                           numPoints = pointPartition->getFiberDimension(*r_iter);
1388:         typename visitor_type::visitor_type nV;
1389:         visitor_type                        cV(*sieve, nV);

1391:         for(int p = 0; p < numPoints; ++p) {
1392:           sieve->cone(points[p], cV);
1393:           if (height) {
1394:             cV.setIsCone(false);
1395:             sieve->support(points[p], cV);
1396:           }
1397:         }
1398:         int i = 0;

1400:         for(typename std::set<typename Mesh::point_type>::const_iterator p_iter = cV.getPoints().begin(); p_iter != cV.getPoints().end(); ++p_iter, ++i) {
1401:           values[i] = *p_iter;
1402:         }
1403:         partition->updatePoint(*r_iter, values);
1404:       }
1405:       delete [] values;
1406:       PETSc::Log::Event("PartitionClosure").end();
1407:     }
1408:     // Create a section mapping points to partitions
1409:     template<typename Section, typename MapSection>
1410:     static void createPartitionMap(const Obj<Section>& partition, const Obj<MapSection>& partitionMap) {
1411:       const typename Section::chart_type& chart = partition->getChart();

1413:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1414:         partitionMap->setFiberDimension(*p_iter, 1);
1415:       }
1416:       partitionMap->allocatePoint();
1417:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1418:         const typename Section::value_type *points = partition->restrictPoint(*p_iter);
1419:         const int                           size   = partition->getFiberDimension(*p_iter);
1420:         const typename Section::point_type  part   = *p_iter;

1422:         for(int i = 0; i < size; ++i) {
1423:           partitionMap->updatePoint(points[i], &part);
1424:         }
1425:       }
1426:     }
1427:   };
1428: #endif

1430:   namespace New {
1431:     template<typename Bundle_, typename Alloc_ = typename Bundle_::alloc_type>
1432:     class Partitioner {
1433:     public:
1434:       typedef Bundle_                          bundle_type;
1435:       typedef Alloc_                           alloc_type;
1436:       typedef typename bundle_type::sieve_type sieve_type;
1437:       typedef typename bundle_type::point_type point_type;
1438:     public:
1441:       // This creates a CSR representation of the adjacency matrix for cells
1442:       // - We allow an exception to contiguous numbering.
1443:       //   If the cell id > numElements, we assign a new number starting at
1444:       //     the top and going downward. I know these might not match up with
1445:       //     the iterator order, but we can fix it later.
1446:       static void buildDualCSR(const Obj<bundle_type>& bundle, const int dim, int **offsets, int **adjacency) {
1447:         ALE_LOG_EVENT_BEGIN;
1448:         typedef typename ALE::New::Completion<bundle_type, point_type, alloc_type> completion;
1449:         const Obj<sieve_type>&                           sieve        = bundle->getSieve();
1450:         const Obj<typename bundle_type::label_sequence>& elements     = bundle->heightStratum(0);
1451:         Obj<sieve_type>                                  overlapSieve = new sieve_type(bundle->comm(), bundle->debug());
1452:         std::map<point_type, point_type>                 newCells;
1453:         int  numElements = elements->size();
1454:         int  newCell     = numElements;
1455:         int *off         = new int[numElements+1];
1456:         int  offset      = 0;
1457:         int *adj;

1459:         completion::scatterSupports(sieve, overlapSieve, bundle->getSendOverlap(), bundle->getRecvOverlap(), bundle);
1460:         if (numElements == 0) {
1461:           *offsets   = NULL;
1462:           *adjacency = NULL;
1463:           ALE_LOG_EVENT_END;
1464:           return;
1465:         }
1466:         if (bundle->depth() == dim) {
1467:           int e = 1;

1469:           off[0] = 0;
1470:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1471:             const Obj<typename sieve_type::traits::coneSequence>& faces  = sieve->cone(*e_iter);
1472:             typename sieve_type::traits::coneSequence::iterator   fBegin = faces->begin();
1473:             typename sieve_type::traits::coneSequence::iterator   fEnd   = faces->end();

1475:             off[e] = off[e-1];
1476:             for(typename sieve_type::traits::coneSequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1477:               if (sieve->support(*f_iter)->size() == 2) {
1478:                 off[e]++;
1479:               } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
1480:                 off[e]++;
1481:               }
1482:             }
1483:             e++;
1484:           }
1485:           adj = new int[off[numElements]];
1486:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1487:             const Obj<typename sieve_type::traits::coneSequence>& faces  = sieve->cone(*e_iter);
1488:             typename sieve_type::traits::coneSequence::iterator   fBegin = faces->begin();
1489:             typename sieve_type::traits::coneSequence::iterator   fEnd   = faces->end();

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

1496:               for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
1497:                 if (*n_iter != *e_iter) adj[offset++] = *n_iter;
1498:               }
1499:               const Obj<typename sieve_type::traits::supportSequence>& oNeighbors = overlapSieve->support(*f_iter);
1500:               typename sieve_type::traits::supportSequence::iterator   onBegin    = oNeighbors->begin();
1501:               typename sieve_type::traits::supportSequence::iterator   onEnd      = oNeighbors->end();

1503:               for(typename sieve_type::traits::supportSequence::iterator n_iter = onBegin; n_iter != onEnd; ++n_iter) {
1504:                 adj[offset++] = *n_iter;
1505:               }
1506:             }
1507:           }
1508:         } else if (bundle->depth() == 1) {
1509:           std::set<point_type> *neighborCells = new std::set<point_type>[numElements];
1510:           int corners      = sieve->cone(*elements->begin())->size();
1511:           int faceVertices = -1;

1513:           if (corners == dim+1) {
1514:             faceVertices = dim;
1515:           } else if ((dim == 2) && (corners == 4)) {
1516:             faceVertices = 2;
1517:           } else if ((dim == 3) && (corners == 8)) {
1518:             faceVertices = 4;
1519:           } else if ((dim == 2) && (corners == 6)) {
1520:             faceVertices = 3;
1521:           } else if ((dim == 2) && (corners == 9)) {
1522:             faceVertices = 3;
1523:           } else if ((dim == 3) && (corners == 10)) {
1524:             faceVertices = 6;
1525:           } else if ((dim == 3) && (corners == 27)) {
1526:             faceVertices = 9;
1527:           } else {
1528:             throw ALE::Exception("Could not determine number of face vertices");
1529:           }
1530:           for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1531:             const Obj<typename sieve_type::traits::coneSequence>& vertices  = sieve->cone(*e_iter);
1532:             typename sieve_type::traits::coneSequence::iterator vEnd = vertices->end();

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

1538:               for(typename sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
1539:                 if (*e_iter == *n_iter) continue;
1540:                 if ((int) sieve->nMeet(*e_iter, *n_iter, 1)->size() == faceVertices) {
1541:                   if ((*e_iter < numElements) && (*n_iter < numElements)) {
1542:                     neighborCells[*e_iter].insert(*n_iter);
1543:                   } else {
1544:                     point_type e = *e_iter, n = *n_iter;

1546:                     if (*e_iter >= numElements) {
1547:                       if (newCells.find(*e_iter) == newCells.end()) newCells[*e_iter] = --newCell;
1548:                       e = newCells[*e_iter];
1549:                     }
1550:                     if (*n_iter >= numElements) {
1551:                       if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
1552:                       n = newCells[*n_iter];
1553:                     }
1554:                     neighborCells[e].insert(n);
1555:                   }
1556:                 }
1557:               }
1558:             }
1559:           }
1560:           off[0] = 0;
1561:           for(int e = 1; e <= numElements; e++) {
1562:             off[e] = neighborCells[e-1].size() + off[e-1];
1563:           }
1564:           adj = new int[off[numElements]];
1565:           for(int e = 0; e < numElements; e++) {
1566:             for(typename std::set<point_type>::iterator n_iter = neighborCells[e].begin(); n_iter != neighborCells[e].end(); ++n_iter) {
1567:               adj[offset++] = *n_iter;
1568:             }
1569:           }
1570:           delete [] neighborCells;
1571:         } else {
1572:           throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
1573:         }
1574:         if (offset != off[numElements]) {
1575:           ostringstream msg;
1576:           msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1577:           throw ALE::Exception(msg.str().c_str());
1578:         }
1579:         //std::cout << "numElements: " << numElements << " newCell: " << newCell << std::endl;
1580:         *offsets   = off;
1581:         *adjacency = adj;
1582:         ALE_LOG_EVENT_END;
1583:       };
1586:       // This creates a CSR representation of the adjacency hypergraph for faces
1587:       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) {
1588:         ALE_LOG_EVENT_BEGIN;
1589:         const Obj<sieve_type>&                           sieve    = bundle->getSieve();
1590:         const Obj<typename bundle_type::label_sequence>& elements = bundle->heightStratum(0);
1591:         int  numElements = elements->size();
1592:         int *off         = new int[numElements+1];
1593:         int  e;

1595:         if (bundle->depth() != dim) {
1596:           throw ALE::Exception("Not yet implemented for non-interpolated meshes");
1597:         }
1598:         off[0] = 0;
1599:         e      = 1;
1600:         for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1601:           off[e] = sieve->cone(*e_iter)->size() + off[e-1];
1602:           e++;
1603:         }
1604:         int *adj    = new int[off[numElements]];
1605:         int  offset = 0;
1606:         for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1607:           const Obj<typename sieve_type::traits::coneSequence>& faces = sieve->cone(*e_iter);
1608:           typename sieve_type::traits::coneSequence::iterator   fEnd  = faces->end();

1610:           for(typename sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
1611:             adj[offset++] = fNumbering->getIndex(*f_iter);
1612:           }
1613:         }
1614:         if (offset != off[numElements]) {
1615:           ostringstream msg;
1616:           msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1617:           throw ALE::Exception(msg.str().c_str());
1618:         }
1619:         *numEdges  = numElements;
1620:         *offsets   = off;
1621:         *adjacency = adj;
1622:         ALE_LOG_EVENT_END;
1623:       };
1624:       template<typename PartitionType>
1625:       static PartitionType *subordinatePartition(const Obj<bundle_type>& bundle, int levels, const Obj<bundle_type>& subBundle, const PartitionType assignment[]) {
1626:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth());
1627:         const Obj<typename bundle_type::label_sequence>& cells      = subBundle->heightStratum(0);
1628:         const Obj<typename bundle_type::numbering_type>& sNumbering = bundle->getFactory()->getLocalNumbering(subBundle, subBundle->depth());
1629:         const int        numCells      = cells->size();
1630:         PartitionType   *subAssignment = new PartitionType[numCells];

1632:         if (levels != 1) {
1633:           throw ALE::Exception("Cannot calculate subordinate partition for any level separation other than 1");
1634:         } else {
1635:           const Obj<typename bundle_type::sieve_type>&   sieve    = bundle->getSieve();
1636:           const Obj<typename bundle_type::sieve_type>&   subSieve = subBundle->getSieve();
1637:           Obj<typename bundle_type::sieve_type::coneSet> tmpSet   = new typename bundle_type::sieve_type::coneSet();
1638:           Obj<typename bundle_type::sieve_type::coneSet> tmpSet2  = new typename bundle_type::sieve_type::coneSet();

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

1643:             Obj<typename bundle_type::sieve_type::supportSet> cell = sieve->nJoin1(cone);
1644:             if (cell->size() != 1) {
1645:               std::cout << "Indeterminate subordinate partition for face " << *c_iter << std::endl;
1646:               for(typename bundle_type::sieve_type::supportSet::iterator s_iter = cell->begin(); s_iter != cell->end(); ++s_iter) {
1647:                 std::cout << "  cell " << *s_iter << std::endl;
1648:               }
1649:               // Could relax this to choosing the first one
1650:               throw ALE::Exception("Indeterminate subordinate partition");
1651:             }
1652:             subAssignment[sNumbering->getIndex(*c_iter)] = assignment[cNumbering->getIndex(*cell->begin())];
1653:             tmpSet->clear();
1654:             tmpSet2->clear();
1655:           }
1656:         }
1657:         return subAssignment;
1658:       }
1659:     };
1660: #ifdef PETSC_HAVE_CHACO
1661:     namespace Chaco {
1662:       template<typename Bundle_>
1663:       class Partitioner {
1664:       public:
1665:         typedef Bundle_                          bundle_type;
1666:         typedef typename bundle_type::sieve_type sieve_type;
1667:         typedef typename bundle_type::point_type point_type;
1668:         typedef short int                        part_type;
1669:       public:
1672:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1673:           part_type *assignment = NULL; /* set number of each vtx (length n) */
1674:           int       *start;             /* start of edge list for each vertex */
1675:           int       *adjacency;         /* = adj -> j; edge list data  */

1677:           ALE_LOG_EVENT_BEGIN;
1678:           ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &start, &adjacency);
1679:           if (bundle->commRank() == 0) {
1680:             /* arguments for Chaco library */
1681:             FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1682:             int nvtxs;                              /* number of vertices in full graph */
1683:             int *vwgts = NULL;                      /* weights for all vertices */
1684:             float *ewgts = NULL;                    /* weights for all edges */
1685:             float *x = NULL, *y = NULL, *z = NULL;  /* coordinates for inertial method */
1686:             char *outassignname = NULL;             /*  name of assignment output file */
1687:             char *outfilename = NULL;               /* output file name */
1688:             int architecture = dim;                 /* 0 => hypercube, d => d-dimensional mesh */
1689:             int ndims_tot = 0;                      /* total number of cube dimensions to divide */
1690:             int mesh_dims[3];                       /* dimensions of mesh of processors */
1691:             double *goal = NULL;                    /* desired set sizes for each set */
1692:             int global_method = 1;                  /* global partitioning algorithm */
1693:             int local_method = 1;                   /* local partitioning algorithm */
1694:             int rqi_flag = 0;                       /* should I use RQI/Symmlq eigensolver? */
1695:             int vmax = 200;                         /* how many vertices to coarsen down to? */
1696:             int ndims = 1;                          /* number of eigenvectors (2^d sets) */
1697:             double eigtol = 0.001;                  /* tolerance on eigenvectors */
1698:             long seed = 123636512;                  /* for random graph mutations */
1699:             float *vCoords[3];

1702:             PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_global_method", &global_method, PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1703:             PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_local_method",  &local_method,  PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1704:             if (global_method == 3) {
1705:               // Inertial Partitioning
1706:               PetscMalloc3(nvtxs,float,&x,nvtxs,float,&y,nvtxs,float,&z);CHKERROR(ierr, "Error in PetscMalloc");
1707:               vCoords[0] = x; vCoords[1] = y; vCoords[2] = z;
1708:               const Obj<typename bundle_type::label_sequence>&    cells       = bundle->heightStratum(0);
1709:               const Obj<typename bundle_type::real_section_type>& coordinates = bundle->getRealSection("coordinates");
1710:               const int corners = bundle->size(coordinates, *(cells->begin()))/dim;
1711:               int       c       = 0;

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

1716:                 for(int d = 0; d < dim; ++d) {
1717:                   vCoords[d][c] = 0.0;
1718:                 }
1719:                 for(int v = 0; v < corners; ++v) {
1720:                   for(int d = 0; d < dim; ++d) {
1721:                     vCoords[d][c] += coords[v*dim+d];
1722:                   }
1723:                 }
1724:                 for(int d = 0; d < dim; ++d) {
1725:                   vCoords[d][c] /= corners;
1726:                 }
1727:               }
1728:             }

1730:             nvtxs = bundle->heightStratum(0)->size();
1731:             mesh_dims[0] = bundle->commSize(); mesh_dims[1] = 1; mesh_dims[2] = 1;
1732:             for(int e = 0; e < start[nvtxs]; e++) {
1733:               adjacency[e]++;
1734:             }
1735:             assignment = new part_type[nvtxs];
1736:             PetscMemzero(assignment, nvtxs * sizeof(part_type));CHKERROR(ierr, "Error in PetscMemzero");

1738:             /* redirect output to buffer: chaco -> msgLog */
1739: #ifdef PETSC_HAVE_UNISTD_H
1740:             char *msgLog;
1741:             int fd_stdout, fd_pipe[2], count;

1743:             fd_stdout = dup(1);
1744:             pipe(fd_pipe);
1745:             close(1);
1746:             dup2(fd_pipe[1], 1);
1747:             msgLog = new char[16284];
1748: #endif

1750:             interface(nvtxs, start, adjacency, vwgts, ewgts, x, y, z,
1751:                              outassignname, outfilename, assignment, architecture, ndims_tot,
1752:                              mesh_dims, goal, global_method, local_method, rqi_flag, vmax, ndims,
1753:                              eigtol, seed);

1755: #ifdef PETSC_HAVE_UNISTD_H
1756:             int SIZE_LOG  = 10000;

1758:             fflush(stdout);
1759:             count = read(fd_pipe[0], msgLog, (SIZE_LOG - 1) * sizeof(char));
1760:             if (count < 0) count = 0;
1761:             msgLog[count] = 0;
1762:             close(1);
1763:             dup2(fd_stdout, 1);
1764:             close(fd_stdout);
1765:             close(fd_pipe[0]);
1766:             close(fd_pipe[1]);
1767:             if (bundle->debug()) {
1768:               std::cout << msgLog << std::endl;
1769:             }
1770:             delete [] msgLog;
1771: #endif
1772:             if (global_method == 3) {
1773:               // Inertial Partitioning
1774:               PetscFree3(x, y, z);CHKERROR(ierr, "Error in PetscFree");
1775:             }
1776:           }
1777:           if (adjacency) delete [] adjacency;
1778:           if (start)     delete [] start;
1779:           ALE_LOG_EVENT_END;
1780:           return assignment;
1781:         };
1782:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1783:           throw ALE::Exception("Chaco cannot partition a mesh by faces");
1784:         };
1785:       };
1786:     };
1787: #endif
1788: #ifdef PETSC_HAVE_PARMETIS
1789:     namespace ParMetis {
1790:       template<typename Bundle_>
1791:       class Partitioner {
1792:       public:
1793:         typedef Bundle_                          bundle_type;
1794:         typedef typename bundle_type::sieve_type sieve_type;
1795:         typedef typename bundle_type::point_type point_type;
1796:         typedef int                              part_type;
1797:       public:
1800:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1801:           int    nvtxs      = 0;    // The number of vertices in full graph
1802:           int   *vtxdist;           // Distribution of vertices across processes
1803:           int   *xadj;              // Start of edge list for each vertex
1804:           int   *adjncy;            // Edge lists for all vertices
1805:           int   *vwgt       = NULL; // Vertex weights
1806:           int   *adjwgt     = NULL; // Edge weights
1807:           int    wgtflag    = 0;    // Indicates which weights are present
1808:           int    numflag    = 0;    // Indicates initial offset (0 or 1)
1809:           int    ncon       = 1;    // The number of weights per vertex
1810:           int    nparts     = bundle->commSize(); // The number of partitions
1811:           float *tpwgts;            // The fraction of vertex weights assigned to each partition
1812:           float *ubvec;             // The balance intolerance for vertex weights
1813:           int    options[5];        // Options
1814:           // Outputs
1815:           int    edgeCut;           // The number of edges cut by the partition
1816:           int   *assignment = NULL; // The vertex partition

1818:           options[0] = 0; // Use all defaults
1819:           vtxdist    = new int[nparts+1];
1820:           vtxdist[0] = 0;
1821:           tpwgts     = new float[ncon*nparts];
1822:           for(int p = 0; p < nparts; ++p) {
1823:             tpwgts[p] = 1.0/nparts;
1824:           }
1825:           ubvec      = new float[ncon];
1826:           ubvec[0]   = 1.05;
1827:           nvtxs      = bundle->heightStratum(0)->size();
1828:           assignment = new part_type[nvtxs];
1829:           MPI_Allgather(&nvtxs, 1, MPI_INT, &vtxdist[1], 1, MPI_INT, bundle->comm());
1830:           for(int p = 2; p <= nparts; ++p) {
1831:             vtxdist[p] += vtxdist[p-1];
1832:           }
1833:           if (bundle->commSize() == 1) {
1834:             PetscMemzero(assignment, nvtxs * sizeof(part_type));
1835:           } else {
1836:             ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &xadj, &adjncy);

1838:             if (bundle->debug() && nvtxs) {
1839:               for(int p = 0; p <= nvtxs; ++p) {
1840:                 std::cout << "["<<bundle->commRank()<<"]xadj["<<p<<"] = " << xadj[p] << std::endl;
1841:               }
1842:               for(int i = 0; i < xadj[nvtxs]; ++i) {
1843:                 std::cout << "["<<bundle->commRank()<<"]adjncy["<<i<<"] = " << adjncy[i] << std::endl;
1844:               }
1845:             }
1846:             if (vtxdist[1] == vtxdist[nparts]) {
1847:               if (bundle->commRank() == 0) {
1848:                 METIS_PartGraphKway(&nvtxs, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &nparts, options, &edgeCut, assignment);
1849:                 if (bundle->debug()) {std::cout << "Metis: edgecut is " << edgeCut << std::endl;}
1850:               }
1851:             } else {
1852:               MPI_Comm comm = bundle->comm();

1854:               ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
1855:               if (bundle->debug()) {std::cout << "ParMetis: edgecut is " << edgeCut << std::endl;}
1856:             }
1857:             if (xadj   != NULL) delete [] xadj;
1858:             if (adjncy != NULL) delete [] adjncy;
1859:           }
1860:           delete [] vtxdist;
1861:           delete [] tpwgts;
1862:           delete [] ubvec;
1863:           return assignment;
1864:         };
1867:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1868: #ifdef PETSC_HAVE_HMETIS
1869:           int   *assignment = NULL; // The vertex partition
1870:           int    nvtxs;      // The number of vertices
1871:           int    nhedges;    // The number of hyperedges
1872:           int   *vwgts;      // The vertex weights
1873:           int   *eptr;       // The offsets of each hyperedge
1874:           int   *eind;       // The vertices in each hyperedge, indexed by eptr
1875:           int   *hewgts;     // The hyperedge weights
1876:           int    nparts;     // The number of partitions
1877:           int    ubfactor;   // The allowed load imbalance (1-50)
1878:           int    options[9]; // Options
1879:           // Outputs
1880:           int    edgeCut;    // The number of edges cut by the partition
1881:           const Obj<ALE::Mesh::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);

1883:           if (topology->commRank() == 0) {
1884:             nvtxs      = bundle->heightStratum(1)->size();
1885:             vwgts      = NULL;
1886:             hewgts     = NULL;
1887:             nparts     = bundle->commSize();
1888:             ubfactor   = 5;
1889:             options[0] = 1;  // Use all defaults
1890:             options[1] = 10; // Number of bisections tested
1891:             options[2] = 1;  // Vertex grouping scheme
1892:             options[3] = 1;  // Objective function
1893:             options[4] = 1;  // V-cycle refinement
1894:             options[5] = 0;
1895:             options[6] = 0;
1896:             options[7] = 1; // Random seed
1897:             options[8] = 24; // Debugging level
1898:             assignment = new part_type[nvtxs];

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

1906:               delete [] eptr;
1907:               delete [] eind;
1908:             }
1909:             if (bundle->debug()) {for (int i = 0; i<nvtxs; i++) printf("[%d] %d\n", PetscGlobalRank, assignment[i]);}
1910:           } else {
1911:             assignment = NULL;
1912:           }
1913:           return assignment;
1914: #else
1915:           throw ALE::Exception("hmetis partitioner is not available.");
1916: #endif
1917:         };
1918:       };
1919:     };
1920: #endif
1921: #ifdef PETSC_HAVE_ZOLTAN
1922:     namespace Zoltan {
1923:       template<typename Bundle_>
1924:       class Partitioner {
1925:       public:
1926:         typedef Bundle_                          bundle_type;
1927:         typedef typename bundle_type::sieve_type sieve_type;
1928:         typedef typename bundle_type::point_type point_type;
1929:         typedef int                              part_type;
1930:       public:
1931:         static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1932:           throw ALE::Exception("Zoltan partition by cells not implemented");
1933:         };
1936:         static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1937:           // Outputs
1938:           float         version;           // The library version
1939:           int           changed;           // Did the partition change?
1940:           int           numGidEntries;     // Number of array entries for a single global ID (1)
1941:           int           numLidEntries;     // Number of array entries for a single local ID (1)
1942:           int           numImport;         // The number of imported points
1943:           ZOLTAN_ID_PTR import_global_ids; // The imported points
1944:           ZOLTAN_ID_PTR import_local_ids;  // The imported points
1945:           int          *import_procs;      // The proc each point was imported from
1946:           int          *import_to_part;    // The partition of each imported point
1947:           int           numExport;         // The number of exported points
1948:           ZOLTAN_ID_PTR export_global_ids; // The exported points
1949:           ZOLTAN_ID_PTR export_local_ids;  // The exported points
1950:           int          *export_procs;      // The proc each point was exported to
1951:           int          *export_to_part;    // The partition assignment of all local points
1952:           int          *assignment;        // The partition assignment of all local points
1953:           const Obj<typename bundle_type::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);

1955:           if (bundle->commSize() == 1) {
1956:             PetscMemzero(assignment, bundle->heightStratum(1)->size() * sizeof(part_type));
1957:           } else {
1958:             if (bundle->commRank() == 0) {
1959:               nvtxs_Zoltan = bundle->heightStratum(1)->size();
1960:               ALE::New::Partitioner<bundle_type>::buildFaceCSR(bundle, dim, fNumbering, &nhedges_Zoltan, &eptr_Zoltan, &eind_Zoltan);
1961:               assignment = new int[nvtxs_Zoltan];
1962:             } else {
1963:               nvtxs_Zoltan   = bundle->heightStratum(1)->size();
1964:               nhedges_Zoltan = 0;
1965:               eptr_Zoltan    = new int[1];
1966:               eind_Zoltan    = new int[1];
1967:               eptr_Zoltan[0] = 0;
1968:               assignment     = NULL;
1969:             }

1971:             int Zoltan_Initialize(0, NULL, &version);
1972:             struct Zoltan_Struct *zz = Zoltan_Create(bundle->comm());
1973:             // General parameters
1974:             Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
1975:             Zoltan_Set_Param(zz, "LB_METHOD", "PHG");
1976:             Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTITION");
1977:             // PHG parameters
1978:             Zoltan_Set_Param(zz, "PHG_OUTPUT_LEVEL", "2");
1979:             Zoltan_Set_Param(zz, "PHG_EDGE_SIZE_THRESHOLD", "1.0"); // Do not throw out dense edges
1980:             // Call backs
1981:             Zoltan_Set_Num_Obj_Fn(zz, getNumVertices_Zoltan, NULL);
1982:             Zoltan_Set_Obj_List_Fn(zz, getLocalElements_Zoltan, NULL);
1983:             Zoltan_Set_HG_Size_CS_Fn(zz, getHgSizes_Zoltan, NULL);
1984:             Zoltan_Set_HG_CS_Fn(zz, getHg_Zoltan, NULL);
1985:             // Debugging
1986:             //Zoltan_Generate_Files(zz, "zoltan.debug", 1, 0, 0, 1); // if using hypergraph callbacks

1988:             Zoltan_LB_Partition(zz, &changed, &numGidEntries, &numLidEntries,
1989:                                        &numImport, &import_global_ids, &import_local_ids, &import_procs, &import_to_part,
1990:                                        &numExport, &export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1991:             for(int v = 0; v < nvtxs_Zoltan; ++v) {
1992:               assignment[v] = export_to_part[v];
1993:             }
1994:             Zoltan_LB_Free_Part(&import_global_ids, &import_local_ids, &import_procs, &import_to_part);
1995:             Zoltan_LB_Free_Part(&export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1996:             Zoltan_Destroy(&zz);

1998:             delete [] eptr_Zoltan;
1999:             delete [] eind_Zoltan;
2000:           }
2001:           if (assignment) {for (int i=0; i<nvtxs_Zoltan; i++) printf("[%d] %d\n",PetscGlobalRank,assignment[i]);}
2002:           return assignment;
2003:         };
2004:       };
2005:     };
2006: #endif
2007:   }
2008: }
2009: #endif