Actual source code: Partitioner.hh
1: #ifndef included_ALE_Partitioner_hh
2: #define included_ALE_Partitioner_hh
4: #ifndef included_ALE_Completion_hh
5: #include <Completion.hh>
6: #endif
8: #ifdef PETSC_HAVE_ZOLTAN
9: #include <zoltan.h>
12: // Inputs
18: int getNumVertices_Zoltan(void *, int *);
20: void getLocalElements_Zoltan(void *, int, int, ZOLTAN_ID_PTR, ZOLTAN_ID_PTR, int, float *, int *);
22: void getHgSizes_Zoltan(void *, int *, int *, int *, int *);
24: void getHg_Zoltan(void *, int, int, int, int, ZOLTAN_ID_PTR, int *, ZOLTAN_ID_PTR, int *);
25: }
27: #endif
29: #ifdef PETSC_HAVE_CHACO
30: #ifdef PETSC_HAVE_UNISTD_H
31: #include <unistd.h>
32: #endif
33: /* Chaco does not have an include file */
36: float *ewgts, float *x, float *y, float *z, char *outassignname,
37: char *outfilename, short *assignment, int architecture, int ndims_tot,
38: int mesh_dims[3], double *goal, int global_method, int local_method,
39: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
42: }
43: #endif
44: #ifdef PETSC_HAVE_PARMETIS
46: #include <parmetis.h>
48: }
49: #endif
50: #ifdef PETSC_HAVE_HMETIS
53: }
54: #endif
56: namespace ALE {
57: #if 1
58: #ifdef PETSC_HAVE_CHACO
59: namespace Chaco {
60: template<typename Alloc_ = malloc_allocator<short int> >
61: class Partitioner {
62: public:
63: typedef short int part_type;
64: typedef Alloc_ alloc_type;
65: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
66: protected:
67: const int logSize;
68: char *msgLog;
69: int fd_stdout, fd_pipe[2];
70: alloc_type _allocator;
71: public:
72: Partitioner(): logSize(10000) {};
73: ~Partitioner() {};
74: protected:
75: // Chaco outputs to stdout. We redirect this to a buffer.
76: // TODO: check error codes for UNIX calls
77: void startStdoutRedirect() {
78: #ifdef PETSC_HAVE_UNISTD_H
79: this->fd_stdout = dup(1);
80: pipe(this->fd_pipe);
81: close(1);
82: dup2(this->fd_pipe[1], 1);
83: #endif
84: };
85: void stopStdoutRedirect() {
86: #ifdef PETSC_HAVE_UNISTD_H
87: int count;
89: fflush(stdout);
90: this->msgLog = new char[this->logSize];
91: count = read(this->fd_pipe[0], this->msgLog, (this->logSize-1)*sizeof(char));
92: if (count < 0) count = 0;
93: this->msgLog[count] = 0;
94: close(1);
95: dup2(this->fd_stdout, 1);
96: close(this->fd_stdout);
97: close(this->fd_pipe[0]);
98: close(this->fd_pipe[1]);
99: //std::cout << this->msgLog << std::endl;
100: delete [] this->msgLog;
101: #endif
102: };
103: public:
104: static bool zeroBase() {return false;}
105: // This method returns the partition section mapping sieve points (here cells) to partitions
106: // start: start of edge list for each vertex
107: // adjacency: adj[start[v]] is edge list data for vertex v
108: // partition: this section is over the partitions and takes points as values
109: // TODO: Read global and local methods from options
110: template<typename Section, typename MeshManager>
111: void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
112: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
113: int nvtxs = numVertices; /* number of vertices in full graph */
114: int *vwgts = NULL; /* weights for all vertices */
115: float *ewgts = NULL; /* weights for all edges */
116: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
117: char *outassignname = NULL; /* name of assignment output file */
118: char *outfilename = NULL; /* output file name */
119: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
120: int ndims_tot = 0; /* total number of cube dimensions to divide */
121: int mesh_dims[3]; /* dimensions of mesh of processors */
122: double *goal = NULL; /* desired set sizes for each set */
123: int global_method = 1; /* global partitioning algorithm */
124: int local_method = 1; /* local partitioning algorithm */
125: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
126: int vmax = 200; /* how many vertices to coarsen down to? */
127: int ndims = 1; /* number of eigenvectors (2^d sets) */
128: double eigtol = 0.001; /* tolerance on eigenvectors */
129: long seed = 123636512; /* for random graph mutations */
130: int maxSize = 0;
132: if (global_method == INERTIAL_METHOD) {manager.createCellCoordinates(nvtxs, &x, &y, &z);}
133: mesh_dims[0] = partition->commSize(); mesh_dims[1] = 1; mesh_dims[2] = 1;
134: part_type *assignment = this->_allocator.allocate(nvtxs);
135: for(int i = 0; i < nvtxs; ++i) {this->_allocator.construct(assignment+i, 0);}
137: this->startStdoutRedirect();
138: interface(nvtxs, (int *) start, (int *) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
139: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
140: vmax, ndims, eigtol, seed);
141: this->stopStdoutRedirect();
143: for(int v = 0; v < nvtxs; ++v) {partition->addFiberDimension(assignment[v], 1);}
144: partition->allocatePoint();
145: for(int p = 0; p < partition->commSize(); ++p) {
146: maxSize = std::max(maxSize, partition->getFiberDimension(p));
147: }
148: typename Section::value_type *values = new typename Section::value_type[maxSize];
150: for(int p = 0; p < partition->commSize(); ++p) {
151: int k = 0;
153: for(int v = 0; v < nvtxs; ++v) {
154: if (assignment[v] == p) values[k++] = manager.getCell(v);
155: }
156: if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
157: partition->updatePoint(p, values);
158: }
159: delete [] values;
161: if (global_method == INERTIAL_METHOD) {manager.destroyCellCoordinates(nvtxs, &x, &y, &z);}
162: for(int i = 0; i < nvtxs; ++i) {this->_allocator.destroy(assignment+i);}
163: this->_allocator.deallocate(assignment, nvtxs);
164: };
165: };
166: }
167: #endif
168: #ifdef PETSC_HAVE_PARMETIS
169: namespace ParMetis {
170: template<typename Alloc_ = malloc_allocator<int> >
171: class Partitioner {
172: public:
173: typedef int part_type;
174: typedef Alloc_ alloc_type;
175: protected:
176: alloc_type _allocator;
177: public:
178: static bool zeroBase() {return true;}
179: // This method returns the partition section mapping sieve points (here cells) to partitions
180: // start: start of edge list for each vertex
181: // adjacency: adj[start[v]] is edge list data for vertex v
182: // partition: this section is over the partitions and takes points as values
183: // TODO: Read parameters from options
184: template<typename Section, typename MeshManager>
185: void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
186: //static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
187: int nvtxs = numVertices; // The number of vertices in full graph
188: int *vtxdist; // Distribution of vertices across processes
189: int *xadj = const_cast<int*>(start); // Start of edge list for each vertex
190: int *adjncy = const_cast<int*>(adjacency); // Edge lists for all vertices
191: int *vwgt = NULL; // Vertex weights
192: int *adjwgt = NULL; // Edge weights
193: int wgtflag = 0; // Indicates which weights are present
194: int numflag = 0; // Indicates initial offset (0 or 1)
195: int ncon = 1; // The number of weights per vertex
196: int nparts = partition->commSize(); // The number of partitions
197: float *tpwgts; // The fraction of vertex weights assigned to each partition
198: float *ubvec; // The balance intolerance for vertex weights
199: int options[5]; // Options
200: int maxSize = 0;
201: // Outputs
202: int edgeCut; // The number of edges cut by the partition
203: part_type *assignment;
205: options[0] = 0; // Use all defaults
206: // Calculate vertex distribution
207: // Not sure this still works in parallel
208: vtxdist = new int[nparts+1];
209: vtxdist[0] = 0;
210: MPI_Allgather(&nvtxs, 1, MPI_INT, &vtxdist[1], 1, MPI_INT, partition->comm());
211: for(int p = 2; p <= nparts; ++p) {
212: vtxdist[p] += vtxdist[p-1];
213: }
214: // Calculate weights
215: tpwgts = new float[ncon*nparts];
216: for(int p = 0; p < nparts; ++p) {
217: tpwgts[p] = 1.0/nparts;
218: }
219: ubvec = new float[ncon];
220: ubvec[0] = 1.05;
222: assignment = this->_allocator.allocate(nvtxs);
223: for(int i = 0; i < nvtxs; ++i) {this->_allocator.construct(assignment+i, 0);}
225: if (partition->commSize() == 1) {
226: PetscMemzero(assignment, nvtxs * sizeof(part_type));
227: } else {
228: if (partition->debug() && nvtxs) {
229: for(int p = 0; p <= nvtxs; ++p) {
230: std::cout << "["<<partition->commRank()<<"]xadj["<<p<<"] = " << xadj[p] << std::endl;
231: }
232: for(int i = 0; i < xadj[nvtxs]; ++i) {
233: std::cout << "["<<partition->commRank()<<"]adjncy["<<i<<"] = " << adjncy[i] << std::endl;
234: }
235: }
236: if (vtxdist[1] == vtxdist[nparts]) {
237: if (partition->commRank() == 0) {
238: METIS_PartGraphKway(&nvtxs, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &nparts, options, &edgeCut, assignment);
239: if (partition->debug()) {std::cout << "Metis: edgecut is " << edgeCut << std::endl;}
240: }
241: } else {
242: MPI_Comm comm = partition->comm();
244: ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
245: if (partition->debug()) {std::cout << "ParMetis: edgecut is " << edgeCut << std::endl;}
246: }
247: }
248: delete [] vtxdist;
249: delete [] tpwgts;
250: delete [] ubvec;
252: for(int v = 0; v < nvtxs; ++v) {partition->addFiberDimension(assignment[v], 1);}
253: partition->allocatePoint();
254: for(int p = 0; p < partition->commSize(); ++p) {
255: maxSize = std::max(maxSize, partition->getFiberDimension(p));
256: }
257: typename Section::value_type *values = new typename Section::value_type[maxSize];
259: for(int p = 0; p < partition->commSize(); ++p) {
260: int k = 0;
262: for(int v = 0; v < nvtxs; ++v) {
263: if (assignment[v] == p) values[k++] = manager.getCell(v);
264: }
265: if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
266: partition->updatePoint(p, values);
267: }
268: delete [] values;
270: for(int i = 0; i < nvtxs; ++i) {this->_allocator.destroy(assignment+i);}
271: this->_allocator.deallocate(assignment, nvtxs);
272: };
273: };
274: };
275: #endif
276: namespace Simple {
277: template<typename Alloc_ = malloc_allocator<short int> >
278: class Partitioner {
279: public:
280: typedef int part_type;
281: typedef Alloc_ alloc_type;
282: protected:
283: alloc_type _allocator;
284: public:
285: Partitioner() {};
286: ~Partitioner() {};
287: public:
288: static bool zeroBase() {return true;}
289: template<typename Section, typename MeshManager>
290: void partition(const int numVertices, const int start[], const int adjacency[], const Obj<Section>& partition, const MeshManager& manager) {
291: const int numProcs = partition->commSize();
292: const int rank = partition->commRank();
293: int maxSize = 0;
295: for(int p = 0; p < numProcs; ++p) {
296: partition->setFiberDimension(p, numVertices/numProcs + ((numVertices % numProcs) > rank));
297: maxSize = std::max(maxSize, partition->getFiberDimension(p));
298: }
299: partition->allocatePoint();
300: typename Section::value_type *values = new typename Section::value_type[maxSize];
302: for(int p = 0; p < partition->commSize(); ++p) {
303: const int start = p*(numVertices/numProcs) + p*((numVertices % numProcs) > p+1);
304: const int end = (p+1)*(numVertices/numProcs) + (p+1)*((numVertices % numProcs) > p+2);
305: int k = 0;
307: for(int v = start; v < end; ++v, ++k) {
308: values[k] = manager.getCell(v);
309: }
310: if (k != partition->getFiberDimension(p)) throw ALE::Exception("Invalid partition");
311: partition->updatePoint(p, values);
312: }
313: delete [] values;
314: };
315: };
316: }
317: #ifdef PETSC_HAVE_CHACO
318: template<typename GraphPartitioner = ALE::Chaco::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
319: #elif defined(PETSC_HAVE_PARMETIS)
320: template<typename GraphPartitioner = ALE::ParMetis::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
321: #else
322: template<typename GraphPartitioner = ALE::Simple::Partitioner<>, typename Alloc_ = malloc_allocator<int> >
323: #endif
324: class Partitioner {
325: public:
326: typedef Alloc_ alloc_type;
327: typedef GraphPartitioner graph_partitioner_type;
328: typedef typename GraphPartitioner::part_type part_type;
329: template<typename Mesh>
330: class MeshManager {
331: public:
332: typedef typename Mesh::point_type point_type;
333: protected:
334: const Obj<Mesh>& mesh;
335: bool simpleCellNumbering;
336: point_type *cells;
337: protected:
338: void createCells(const int height) {
339: const Obj<typename Mesh::label_sequence>& mcells = mesh->heightStratum(height);
340: const typename Mesh::label_sequence::iterator cEnd = mcells->end();
341: const int numCells = mcells->size();
342: int c = 0;
344: this->cells = NULL;
345: this->simpleCellNumbering = true;
346: for(typename Mesh::label_sequence::iterator c_iter = mcells->begin(); c_iter != cEnd; ++c_iter, ++c) {
347: if (*c_iter != c) {
348: this->simpleCellNumbering = false;
349: break;
350: }
351: }
352: if (!this->simpleCellNumbering) {
353: this->cells = new point_type[numCells];
354: c = 0;
355: for(typename Mesh::label_sequence::iterator c_iter = mcells->begin(); c_iter != cEnd; ++c_iter, ++c) {
356: this->cells[c] = *c_iter;
357: }
358: }
359: };
360: public:
361: MeshManager(const Obj<Mesh>& mesh, const int height = 0): mesh(mesh) {
362: this->createCells(height);
363: };
364: ~MeshManager() {
365: if (this->cells) {delete [] this->cells;}
366: };
367: public:
368: template<typename Float>
369: void createCellCoordinates(const int numVertices, Float *X[], Float *Y[], Float *Z[]) const {
370: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
371: const int dim = mesh->getDimension();
372: typedef typename alloc_type::template rebind<Float>::other float_alloc_type;
373: Float *x = float_alloc_type().allocate(numVertices*3);
374: for(int i = 0; i < numVertices*3; ++i) {float_alloc_type().construct(x+i, 0.0);}
375: Float *y = x+numVertices;
376: Float *z = y+numVertices;
377: Float *vCoords[3];
379: vCoords[0] = x; vCoords[1] = y; vCoords[2] = z;
380: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
381: const int corners = mesh->size(coordinates, *(cells->begin()))/dim;
382: int c = 0;
384: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter !=cells->end(); ++c_iter, ++c) {
385: const double *coords = mesh->restrictClosure(coordinates, *c_iter);
387: for(int d = 0; d < dim; ++d) {
388: vCoords[d][c] = 0.0;
389: }
390: for(int v = 0; v < corners; ++v) {
391: for(int d = 0; d < dim; ++d) {
392: vCoords[d][c] += coords[v*dim+d];
393: }
394: }
395: for(int d = 0; d < dim; ++d) {
396: vCoords[d][c] /= corners;
397: }
398: }
399: *X = x;
400: *Y = y;
401: *Z = z;
402: };
403: template<typename Float>
404: void destroyCellCoordinates(const int numVertices, Float *X[], Float *Y[], Float *Z[]) const {
405: typedef typename alloc_type::template rebind<Float>::other float_alloc_type;
406: Float *x = *X;
408: for(int i = 0; i < numVertices*3; ++i) {float_alloc_type().destroy(x+i);}
409: float_alloc_type().deallocate(x, numVertices*3);
410: };
411: point_type getCell(const int cellNumber) const {
412: if (this->simpleCellNumbering) {
413: return cellNumber;
414: }
415: return this->cells[cellNumber];
416: };
417: };
418: template<typename Sieve>
419: class OffsetVisitor {
420: const Sieve& sieve;
421: const Sieve& overlapSieve;
422: int *offsets;
423: public:
424: OffsetVisitor(const Sieve& s, const Sieve& ovS, int off[]) : sieve(s), overlapSieve(ovS), offsets(off) {};
425: void visitPoint(const typename Sieve::point_type& point) {};
426: void visitArrow(const typename Sieve::arrow_type& arrow) {
427: const typename Sieve::point_type cell = arrow.target;
428: const typename Sieve::point_type face = arrow.source;
429: const int size = this->sieve.getSupportSize(face);
430: const int ovSize = this->overlapSieve.getSupportSize(face);
432: if (size == 2) {
433: offsets[cell+1]++;
434: } else if ((size == 1) && (ovSize == 1)) {
435: offsets[cell+1]++;
436: }
437: };
438: };
439: template<typename Sieve>
440: class AdjVisitor {
441: protected:
442: typename Sieve::point_type cell;
443: int *adjacency;
444: const int cellOffset;
445: int offset;
446: public:
447: AdjVisitor(int adj[], const bool zeroBase) : adjacency(adj), cellOffset(zeroBase ? 0 : 1), offset(0) {};
448: void visitPoint(const typename Sieve::point_type& point) {};
449: void visitArrow(const typename Sieve::arrow_type& arrow) {
450: const int neighbor = arrow.target;
452: if (neighbor != this->cell) {
453: //std::cout << "Adding dual edge from " << cell << " to " << neighbor << std::endl;
454: this->adjacency[this->offset++] = neighbor + this->cellOffset;
455: }
456: };
457: public:
458: void setCell(const typename Sieve::point_type cell) {this->cell = cell;};
459: int getOffset() {return this->offset;}
460: };
461: template<typename Sieve>
462: class MeetVisitor {
463: public:
464: typedef std::set<typename Sieve::point_type> neighbors_type;
465: protected:
466: const Sieve& sieve;
467: const int numCells;
468: const int faceVertices;
469: typename Sieve::point_type newCell;
470: neighbors_type *neighborCells;
471: typename ISieveVisitor::PointRetriever<Sieve> *pR;
472: typename Sieve::point_type cell;
473: std::map<typename Sieve::point_type, typename Sieve::point_type> newCells;
474: public:
475: MeetVisitor(const Sieve& s, const int n, const int fV) : sieve(s), numCells(n), faceVertices(fV), newCell(n) {
476: this->neighborCells = new std::set<typename Sieve::point_type>[numCells];
477: this->pR = new typename ISieveVisitor::PointRetriever<Sieve>(this->sieve.getMaxConeSize());
478: };
479: ~MeetVisitor() {delete [] this->neighborCells; delete this->pR;};
480: void visitArrow(const typename Sieve::arrow_type& arrow) {};
481: void visitPoint(const typename Sieve::point_type& point) {
482: const typename Sieve::point_type& neighbor = point;
484: if (this->cell == neighbor) return;
485: this->pR->clear();
486: this->sieve.meet(this->cell, neighbor, *this->pR);
487: if (this->pR->getSize() == (size_t) this->faceVertices) {
488: if ((this->cell < numCells) && (neighbor < numCells)) {
489: this->neighborCells[this->cell].insert(neighbor);
490: } else {
491: typename Sieve::point_type e = this->cell, n = neighbor;
493: if (this->cell >= numCells) {
494: if (this->newCells.find(cell) == this->newCells.end()) this->newCells[cell] = --newCell;
495: e = this->newCells[cell];
496: }
497: if (neighbor >= numCells) {
498: if (this->newCells.find(neighbor) == this->newCells.end()) this->newCells[neighbor] = --newCell;
499: n = this->newCells[neighbor];
500: }
501: this->neighborCells[e].insert(n);
502: }
503: }
504: };
505: public:
506: void setCell(const typename Sieve::point_type& c) {this->cell = c;};
507: const neighbors_type *getNeighbors() {return this->neighborCells;};
508: };
509: public: // Creating overlaps
510: // Create a partition point overlap for distribution
511: // This is the default overlap which comes from distributing a serial mesh on process 0
512: template<typename SendOverlap, typename RecvOverlap>
513: static void createDistributionPartOverlap(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
514: const int rank = sendOverlap->commRank();
515: const int size = sendOverlap->commSize();
517: if (rank == 0) {
518: for(int p = 1; p < size; p++) {
519: // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
520: sendOverlap->addCone(p, p, p);
521: }
522: }
523: if (rank != 0) {
524: // The arrow is from remote partition point rank (color) on rank 0 (source) to local partition point rank (target)
525: recvOverlap->addCone(0, rank, rank);
526: }
527: };
528: // Create a mesh point overlap for distribution
529: // A local numbering is created for the remote points
530: // This is the default overlap which comes from distributing a serial mesh on process 0
531: template<typename Section, typename RecvPartOverlap, typename Renumbering, typename SendOverlap, typename RecvOverlap>
532: static void createDistributionMeshOverlap(const Obj<Section>& partition, const Obj<RecvPartOverlap>& recvPartOverlap, Renumbering& renumbering, const Obj<Section>& overlapPartition, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
533: const typename Section::chart_type& chart = partition->getChart();
535: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
536: if (*p_iter == sendOverlap->commRank()) continue;
537: const typename Section::value_type *points = partition->restrictPoint(*p_iter);
538: const int numPoints = partition->getFiberDimension(*p_iter);
539:
540: for(int i = 0; i < numPoints; ++i) {
541: // Notice here that we do not know the local renumbering (but we do not use it)
542: sendOverlap->addArrow(points[i], *p_iter, points[i]);
543: }
544: }
545: if (sendOverlap->debug()) {sendOverlap->view("Send mesh overlap");}
546: const Obj<typename RecvPartOverlap::traits::baseSequence> rPoints = recvPartOverlap->base();
548: for(typename RecvPartOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
549: const Obj<typename RecvPartOverlap::coneSequence>& ranks = recvPartOverlap->cone(*p_iter);
550: //const typename Section::point_type& localPartPoint = *p_iter;
551: const typename Section::point_type rank = *ranks->begin();
552: const typename Section::point_type& remotePartPoint = ranks->begin().color();
553: const typename Section::value_type *points = overlapPartition->restrictPoint(remotePartPoint);
554: const int numPoints = overlapPartition->getFiberDimension(remotePartPoint);
556: for(int i = 0; i < numPoints; ++i) {
557: recvOverlap->addArrow(rank, renumbering[points[i]], points[i]);
558: }
559: }
560: if (recvOverlap->debug()) {recvOverlap->view("Receive mesh overlap");}
561: };
562: // Create a partition point overlap from a partition
563: // The intention is to create an overlap which enables exchange of redistribution information
564: template<typename Section, typename SendOverlap, typename RecvOverlap>
565: static void createPartitionPartOverlap(const Obj<Section>& partition, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
566: const typename Section::chart_type& chart = partition->getChart();
567: const int rank = partition->commRank();
568: const int size = partition->commSize();
569: int *adj = new int[size];
570: int *recvCounts = new int[size];
571: int numNeighbors;
573: for(int p = 0; p < size; ++p) {
574: adj[p] = 0;
575: recvCounts[p] = 1;
576: }
577: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart->end(); ++p_iter) {
578: const typename Section::value_type& p = partition->restrictPoint(*p_iter)[0];
579: // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
580: sendOverlap->addCone(p, p, p);
581: adj[p] = 1;
582: }
583: MPI_Reduce_Scatter(adj, &numNeighbors, recvCounts, size, MPI_INT, MPI_SUM, partition->comm());
584: MPI_Request *recvRequests = new MPI_Request[numNeighbors];
585: int dummy = 0;
587: // TODO: Get a unique tag
588: for(int n = 0; n < numNeighbors; ++n) {
589: MPI_Irecv(&dummy, 1, MPI_INT, MPI_ANY_SOURCE, 1, partition->comm(), &recvRequests[n]);
590: }
591: const Obj<typename SendOverlap::traits::baseSequence> ranks = sendOverlap->base();
592: const typename SendOverlap::traits::baseSequence::iterator rEnd = ranks->end();
593: MPI_Request *sendRequests = new MPI_Request[ranks->size()];
594: int s = 0;
596: for(typename SendOverlap::traits::baseSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter, ++s) {
597: MPI_Isend(&dummy, 1, MPI_INT, *r_iter, 1, partition->comm(), &sendRequests[s]);
598: }
599: for(int n = 0; n < numNeighbors; ++n) {
600: MPI_Status status;
601: int idx;
603: MPI_Waitany(numNeighbors, recvRequests, &idx, &status);
604: // The arrow is from remote partition point rank (color) on rank p (source) to local partition point rank (target)
605: recvOverlap->addCone(status.MPI_SOURCE, rank, rank);
606: }
607: MPI_Waitall(ranks->size(), sendRequests, MPI_STATUSES_IGNORE);
608: delete [] sendRequests;
609: delete [] recvRequests;
610: delete [] adj;
611: delete [] recvCounts;
612: };
613: public: // Building CSR meshes
614: // This produces the dual graph (each cell is a vertex and each face is an edge)
615: // numbering: A contiguous numbering of the cells (not yet included)
616: // numVertices: The number of vertices in the graph (cells in the mesh)
617: // adjacency: The vertices adjacent to each vertex (cells adjacent to each mesh cell)
618: // - We allow an exception to contiguous numbering.
619: // If the cell id > numElements, we assign a new number starting at
620: // the top and going downward. I know these might not match up with
621: // the iterator order, but we can fix it later.
622: template<typename Mesh>
623: static void buildDualCSR(const Obj<Mesh>& mesh, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
624: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
625: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
626: const typename Mesh::label_sequence::iterator cEnd = cells->end();
627: const int numCells = cells->size();
628: int newCell = numCells;
629: Obj<typename Mesh::sieve_type> overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
630: int offset = 0;
631: const int cellOffset = zeroBase ? 0 : 1;
632: const int dim = mesh->getDimension();
633: std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;
635: // TODO: This is necessary for parallel partitioning
636: //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
637: if (numCells == 0) {
638: *numVertices = 0;
639: *offsets = NULL;
640: *adjacency = NULL;
641: return;
642: }
643: int *off = alloc_type().allocate(numCells+1);
644: int *adj;
645: for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
646: if (mesh->depth() == dim) {
647: int c = 1;
649: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
650: const Obj<typename Mesh::sieve_type::traits::coneSequence>& faces = sieve->cone(*c_iter);
651: const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd = faces->end();
653: off[c] = off[c-1];
654: for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
655: if (sieve->support(*f_iter)->size() == 2) {
656: off[c]++;
657: } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
658: off[c]++;
659: }
660: }
661: }
662: adj = alloc_type().allocate(off[numCells]);
663: for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
664: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
665: const Obj<typename Mesh::sieve_type::traits::coneSequence>& faces = sieve->cone(*c_iter);
666: const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd = faces->end();
668: for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
669: const Obj<typename Mesh::sieve_type::traits::supportSequence>& neighbors = sieve->support(*f_iter);
670: const typename Mesh::sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
672: for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
673: if (*n_iter != *c_iter) adj[offset++] = *n_iter + cellOffset;
674: }
675: const Obj<typename Mesh::sieve_type::traits::supportSequence>& oNeighbors = overlapSieve->support(*f_iter);
676: const typename Mesh::sieve_type::traits::supportSequence::iterator onEnd = oNeighbors->end();
678: for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = oNeighbors->begin(); n_iter != onEnd; ++n_iter) {
679: adj[offset++] = *n_iter + cellOffset;
680: }
681: }
682: }
683: } else if (mesh->depth() == 1) {
684: std::set<typename Mesh::point_type> *neighborCells = new std::set<typename Mesh::point_type>[numCells];
685: const int corners = sieve->cone(*cells->begin())->size();
686: int faceVertices;
688: if (corners == dim+1) {
689: faceVertices = dim;
690: } else if ((dim == 2) && (corners == 4)) {
691: faceVertices = 2;
692: } else if ((dim == 3) && (corners == 8)) {
693: faceVertices = 4;
694: } else {
695: throw ALE::Exception("Could not determine number of face vertices");
696: }
697: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
698: const Obj<typename Mesh::sieve_type::traits::coneSequence>& vertices = sieve->cone(*c_iter);
699: const typename Mesh::sieve_type::traits::coneSequence::iterator vEnd = vertices->end();
701: for(typename Mesh::sieve_type::traits::coneSequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
702: const Obj<typename Mesh::sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
703: const typename Mesh::sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
705: for(typename Mesh::sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
706: if (*c_iter == *n_iter) continue;
707: if ((int) sieve->nMeet(*c_iter, *n_iter, 1)->size() == faceVertices) {
708: if ((*c_iter < numCells) && (*n_iter < numCells)) {
709: neighborCells[*c_iter].insert(*n_iter);
710: } else {
711: typename Mesh::point_type e = *c_iter, n = *n_iter;
713: if (*c_iter >= numCells) {
714: if (newCells.find(*c_iter) == newCells.end()) newCells[*c_iter] = --newCell;
715: e = newCells[*c_iter];
716: }
717: if (*n_iter >= numCells) {
718: if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
719: n = newCells[*n_iter];
720: }
721: neighborCells[e].insert(n);
722: }
723: }
724: }
725: }
726: }
727: off[0] = 0;
728: for(int c = 1; c <= numCells; c++) {
729: off[c] = neighborCells[c-1].size() + off[c-1];
730: }
731: adj = alloc_type().allocate(off[numCells]);
732: for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
733: for(int c = 0; c < numCells; c++) {
734: for(typename std::set<typename Mesh::point_type>::iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
735: adj[offset++] = *n_iter + cellOffset;
736: }
737: }
738: delete [] neighborCells;
739: } else {
740: throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
741: }
742: if (offset != off[numCells]) {
743: ostringstream msg;
744: msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
745: throw ALE::Exception(msg.str().c_str());
746: }
747: *numVertices = numCells;
748: *offsets = off;
749: *adjacency = adj;
750: };
751: template<typename Mesh>
752: static void buildDualCSRV(const Obj<Mesh>& mesh, int *numVertices, int **offsets, int **adjacency, const bool zeroBase = true) {
753: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
754: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
755: const typename Mesh::label_sequence::iterator cEnd = cells->end();
756: const int numCells = cells->size();
757: Obj<typename Mesh::sieve_type> overlapSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
758: int offset = 0;
759: const int cellOffset = zeroBase ? 0 : 1;
760: const int dim = mesh->getDimension();
761: std::map<typename Mesh::point_type, typename Mesh::point_type> newCells;
763: // TODO: This is necessary for parallel partitioning
764: //completion::scatterSupports(sieve, overlapSieve, mesh->getSendOverlap(), mesh->getRecvOverlap(), mesh);
765: overlapSieve->setChart(sieve->getChart());
766: overlapSieve->allocate();
767: if (numCells == 0) {
768: *numVertices = 0;
769: *offsets = NULL;
770: *adjacency = NULL;
771: return;
772: }
773: int *off = alloc_type().allocate(numCells+1);
774: int *adj;
775: for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
776: if (mesh->depth() == dim) {
777: OffsetVisitor<typename Mesh::sieve_type> oV(*sieve, *overlapSieve, off);
779: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
780: sieve->cone(*c_iter, oV);
781: }
782: for(int p = 1; p <= numCells; ++p) {
783: off[p] = off[p] + off[p-1];
784: }
785: adj = alloc_type().allocate(off[numCells]);
786: for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
787: AdjVisitor<typename Mesh::sieve_type> aV(adj, zeroBase);
788: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > sV(*sieve, aV);
789: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > ovSV(*overlapSieve, aV);
791: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
792: aV.setCell(*c_iter);
793: sieve->cone(*c_iter, sV);
794: sieve->cone(*c_iter, ovSV);
795: }
796: offset = aV.getOffset();
797: } else if (mesh->depth() == 1) {
798: typedef MeetVisitor<typename Mesh::sieve_type> mv_type;
799: typedef typename ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, mv_type> sv_type;
800: const int corners = sieve->getConeSize(*cells->begin());
801: int faceVertices;
803: if (corners == dim+1) {
804: faceVertices = dim;
805: } else if ((dim == 2) && (corners == 4)) {
806: faceVertices = 2;
807: } else if ((dim == 3) && (corners == 8)) {
808: faceVertices = 4;
809: } else {
810: throw ALE::Exception("Could not determine number of face vertices");
811: }
812: mv_type mV(*sieve, numCells, faceVertices);
813: sv_type sV(*sieve, mV);
815: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
816: mV.setCell(*c_iter);
817: sieve->cone(*c_iter, sV);
818: }
819: const typename mv_type::neighbors_type *neighborCells = mV.getNeighbors();
821: off[0] = 0;
822: for(int c = 1; c <= numCells; c++) {
823: off[c] = neighborCells[c-1].size() + off[c-1];
824: }
825: adj = alloc_type().allocate(off[numCells]);
826: for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
827: for(int c = 0; c < numCells; c++) {
828: for(typename mv_type::neighbors_type::const_iterator n_iter = neighborCells[c].begin(); n_iter != neighborCells[c].end(); ++n_iter) {
829: //std::cout << "Adding dual edge from " << c << " to " << *n_iter << std::endl;
830: adj[offset++] = *n_iter + cellOffset;
831: }
832: }
833: } else {
834: throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
835: }
836: if (offset != off[numCells]) {
837: ostringstream msg;
838: msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
839: throw ALE::Exception(msg.str().c_str());
840: }
841: *numVertices = numCells;
842: *offsets = off;
843: *adjacency = adj;
844: };
845: // This produces a hypergraph (each face is a vertex and each cell is a hyperedge)
846: // numbering: A contiguous numbering of the faces
847: // numEdges: The number of edges in the hypergraph
848: // adjacency: The vertices in each edge
849: template<typename Mesh>
850: static void buildFaceDualCSR(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
851: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
852: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
853: const typename Mesh::label_sequence::iterator cEnd = cells->end();
854: const int faceOffset = zeroBase ? 0 : 1;
855: int numCells = cells->size();
856: int c = 1;
858: if (mesh->depth() != mesh->getDimension()) {throw ALE::Exception("Not yet implemented for non-interpolated meshes");}
859: int *off = alloc_type().allocate(numCells+1);
860: for(int i = 0; i < numCells+1; ++i) {alloc_type().construct(off+i, 0);}
861: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
862: off[c] = sieve->cone(*c_iter)->size() + off[c-1];
863: }
864: int *adj = alloc_type().allocate(off[numCells]);
865: for(int i = 0; i < off[numCells]; ++i) {alloc_type().construct(adj+i, 0);}
866: int offset = 0;
867: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
868: const Obj<typename Mesh::sieve_type::traits::coneSequence>& faces = sieve->cone(*c_iter);
869: const typename Mesh::sieve_type::traits::coneSequence::iterator fEnd = faces->end();
871: for(typename Mesh::sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter, ++offset) {
872: adj[offset] = numbering->getIndex(*f_iter) + faceOffset;
873: }
874: }
875: if (offset != off[numCells]) {
876: ostringstream msg;
877: msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numCells];
878: throw ALE::Exception(msg.str().c_str());
879: }
880: *numEdges = numCells;
881: *offsets = off;
882: *adjacency = adj;
883: };
884: template<typename Mesh>
885: static void buildFaceDualCSRV(const Obj<Mesh>& mesh, const Obj<typename Mesh::numbering_type>& numbering, int *numEdges, int **offsets, int **adjacency, const bool zeroBase = true) {
886: throw ALE::Exception("Not implemented");
887: };
888: static void destroyCSR(int numPoints, int *offsets, int *adjacency) {
889: if (adjacency) {
890: for(int i = 0; i < offsets[numPoints]; ++i) {alloc_type().destroy(adjacency+i);}
891: alloc_type().deallocate(adjacency, offsets[numPoints]);
892: }
893: if (offsets) {
894: for(int i = 0; i < numPoints+1; ++i) {alloc_type().destroy(offsets+i);}
895: alloc_type().deallocate(offsets, numPoints+1);
896: }
897: };
898: template<typename OldSection, typename Partition, typename Renumbering, typename NewSection>
899: static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<NewSection>& newSection) {
900: const typename Partition::value_type *points = partition->restrictPoint(oldSection->commRank());
901: const int numPoints = partition->getFiberDimension(oldSection->commRank());
903: for(int p = 0; p < numPoints; ++p) {
904: if (oldSection->hasPoint(points[p])) {
905: newSection->setFiberDimension(renumbering[points[p]], oldSection->getFiberDimension(points[p]));
906: }
907: }
908: if (numPoints) {newSection->allocatePoint();}
909: for(int p = 0; p < numPoints; ++p) {
910: if (oldSection->hasPoint(points[p])) {
911: newSection->updatePointAll(renumbering[points[p]], oldSection->restrictPoint(points[p]));
912: }
913: }
914: };
915: // Specialize to ArrowSection
916: template<typename OldSection, typename Partition, typename Renumbering>
917: static void createLocalSection(const Obj<OldSection>& oldSection, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<UniformSection<MinimalArrow<int,int>,int> >& newSection) {
918: typedef UniformSection<MinimalArrow<int,int>,int> NewSection;
919: const typename Partition::value_type *points = partition->restrictPoint(oldSection->commRank());
920: const int numPoints = partition->getFiberDimension(oldSection->commRank());
921: const typename OldSection::chart_type& oldChart = oldSection->getChart();
922: std::set<typename Partition::value_type> myPoints;
924: for(int p = 0; p < numPoints; ++p) {
925: myPoints.insert(points[p]);
926: }
927: for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
928: if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
929: newSection->setFiberDimension(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), oldSection->getFiberDimension(*c_iter));
930: }
931: }
932: if (oldChart.size()) {newSection->allocatePoint();}
933: for(typename OldSection::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
934: if (myPoints.count(c_iter->source) && myPoints.count(c_iter->target)) {
935: const typename OldSection::value_type *values = oldSection->restrictPoint(*c_iter);
937: newSection->updatePointAll(typename OldSection::point_type(renumbering[c_iter->source], renumbering[c_iter->target]), values);
938: }
939: }
940: };
941: template<typename Sifter, typename Section, typename Renumbering>
942: static void createLocalSifter(const Obj<Sifter>& sifter, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sifter>& localSifter) {
943: const typename Section::value_type *points = partition->restrictPoint(sifter->commRank());
944: const int numPoints = partition->getFiberDimension(sifter->commRank());
946: for(int p = 0; p < numPoints; ++p) {
947: const Obj<typename Sifter::traits::coneSequence>& cone = sifter->cone(points[p]);
948: const typename Sifter::traits::coneSequence::iterator cEnd = cone->end();
950: for(typename Sifter::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
951: localSifter->addArrow(*c_iter, renumbering[points[p]]);
952: }
953: }
954: };
955: template<typename Sieve, typename Section, typename Renumbering>
956: static void createLocalSieve(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
957: const typename Section::value_type *points = partition->restrictPoint(sieve->commRank());
958: const int numPoints = partition->getFiberDimension(sieve->commRank());
960: for(int p = 0; p < numPoints; ++p) {
961: Obj<typename Sieve::coneSet> current = new typename Sieve::coneSet();
962: Obj<typename Sieve::coneSet> next = new typename Sieve::coneSet();
963: Obj<typename Sieve::coneSet> tmp;
965: current->insert(points[p]);
966: while(current->size()) {
967: for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
968: const Obj<typename Sieve::traits::coneSequence>& cone = sieve->cone(*p_iter);
969:
970: for(typename Sieve::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
971: localSieve->addArrow(renumbering[*c_iter], renumbering[*p_iter], c_iter.color());
972: next->insert(*c_iter);
973: }
974: }
975: tmp = current; current = next; next = tmp;
976: next->clear();
977: }
978: if (height) {
979: current->insert(points[p]);
980: while(current->size()) {
981: for(typename Sieve::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
982: const Obj<typename Sieve::traits::supportSequence>& support = sieve->support(*p_iter);
983:
984: for(typename Sieve::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
985: localSieve->addArrow(renumbering[*p_iter], renumbering[*s_iter], s_iter.color());
986: next->insert(*s_iter);
987: }
988: }
989: tmp = current; current = next; next = tmp;
990: next->clear();
991: }
992: }
993: }
994: };
995: template<typename Mesh, typename Section, typename Renumbering>
996: static void createLocalMesh(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
997: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
998: const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();
1000: createLocalSieve(sieve, partition, renumbering, localSieve, height);
1001: };
1002: template<typename Sieve, typename Section, typename Renumbering>
1003: static void sizeLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1004: typedef std::set<typename Sieve::point_type> pointSet;
1005: const typename Section::value_type *points = partition->restrictPoint(sieve->commRank());
1006: const int numPoints = partition->getFiberDimension(sieve->commRank());
1007: int maxSize = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1008: const pointSet pSet(points, &points[numPoints]);
1009: ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> fV(pSet, renumbering, maxSize);
1011: for(int p = 0; p < numPoints; ++p) {
1012: sieve->cone(points[p], fV);
1013: localSieve->setConeSize(renumbering[points[p]], fV.getSize());
1014: fV.clear();
1015: sieve->support(points[p], fV);
1016: localSieve->setSupportSize(renumbering[points[p]], fV.getSize());
1017: fV.clear();
1018: }
1019: };
1020: template<typename Mesh, typename Section, typename Renumbering>
1021: static void sizeLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1022: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1023: const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();
1025: sizeLocalSieveV(sieve, partition, renumbering, localSieve, height);
1026: };
1027: template<typename Sieve, typename Section, typename Renumbering>
1028: static void createLocalLabelV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1029: typedef std::set<typename Sieve::point_type> pointSet;
1030: typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1031: const typename Section::value_type *points = partition->restrictPoint(sieve->commRank());
1032: const int numPoints = partition->getFiberDimension(sieve->commRank());
1033: int maxSize = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1034: typename Sieve::point_type *oPoints = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1035: int *oOrients = new int[std::max(1, sieve->getMaxConeSize())];
1036: const pointSet pSet(points, &points[numPoints]);
1037: visitor_type fV(pSet, renumbering, maxSize);
1039: for(int p = 0; p < numPoints; ++p) {
1040: fV.useRenumbering(false);
1041: sieve->orientedCone(points[p], fV);
1042: const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1043: const int n = fV.getOrientedSize();
1044: for(int i = 0; i < n; ++i) {
1045: oPoints[i] = q[i].first;
1046: oOrients[i] = q[i].second;
1047: }
1048: localSieve->setCone(oPoints, renumbering[points[p]]);
1049: localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1050: fV.clear();
1051: fV.useRenumbering(true);
1052: sieve->support(points[p], fV);
1053: if (fV.getSize()) {localSieve->setSupport(points[p], fV.getPoints());}
1054: fV.clear();
1055: }
1056: delete [] oPoints;
1057: delete [] oOrients;
1058: };
1059: template<typename Sieve, typename Section, typename Renumbering>
1060: static void createLocalSieveV(const Obj<Sieve>& sieve, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Sieve>& localSieve, const int height = 0) {
1061: typedef std::set<typename Sieve::point_type> pointSet;
1062: typedef ISieveVisitor::FilteredPointRetriever<Sieve,pointSet,Renumbering> visitor_type;
1063: const typename Section::value_type *points = partition->restrictPoint(sieve->commRank());
1064: const int numPoints = partition->getFiberDimension(sieve->commRank());
1065: int maxSize = std::max(0, std::max(sieve->getMaxConeSize(), sieve->getMaxSupportSize()));
1066: typename Sieve::point_type *oPoints = new typename Sieve::point_type[std::max(1, sieve->getMaxConeSize())];
1067: int *oOrients = new int[std::max(1, sieve->getMaxConeSize())];
1068: const pointSet pSet(points, &points[numPoints]);
1069: visitor_type fV(pSet, renumbering, maxSize);
1071: for(int p = 0; p < numPoints; ++p) {
1072: ///sieve->cone(points[p], fV);
1073: ///localSiaeve->setCone(fV.getPoints(), renumbering[points[p]]);
1074: sieve->orientedCone(points[p], fV);
1075: const typename visitor_type::oriented_point_type *q = fV.getOrientedPoints();
1076: const int n = fV.getOrientedSize();
1077: for(int i = 0; i < n; ++i) {
1078: oPoints[i] = q[i].first;
1079: oOrients[i] = q[i].second;
1080: }
1081: localSieve->setCone(oPoints, renumbering[points[p]]);
1082: localSieve->setConeOrientation(oOrients, renumbering[points[p]]);
1083: fV.clear();
1084: sieve->support(points[p], fV);
1085: localSieve->setSupport(renumbering[points[p]], fV.getPoints());
1086: fV.clear();
1087: }
1088: delete [] oPoints;
1089: delete [] oOrients;
1090: };
1091: template<typename Mesh, typename Section, typename Renumbering>
1092: static void createLocalMeshV(const Obj<Mesh>& mesh, const Obj<Section>& partition, Renumbering& renumbering, const Obj<Mesh>& localMesh, const int height = 0) {
1093: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1094: const Obj<typename Mesh::sieve_type>& localSieve = localMesh->getSieve();
1096: createLocalSieveV(sieve, partition, renumbering, localSieve, height);
1097: };
1098: public: // Partitioning
1099: // partition: Should be properly allocated on input
1100: // height: Height of the point set to uniquely partition
1101: // TODO: Could rebind assignment section to the type of the output
1102: template<typename Mesh, typename Section>
1103: static void createPartition(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1104: MeshManager<Mesh> manager(mesh, height);
1105: int *start = NULL;
1106: int *adjacency = NULL;
1108: if (height == 0) {
1109: int numVertices;
1111: buildDualCSR(mesh, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1112: GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1113: destroyCSR(numVertices, start, adjacency);
1114: } else if (height == 1) {
1115: int numEdges;
1117: buildFaceDualCSR(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1118: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1119: destroyCSR(numEdges, start, adjacency);
1120: } else {
1121: throw ALE::Exception("Invalid partition height");
1122: }
1123: };
1124: template<typename Mesh, typename Section>
1125: static void createPartitionV(const Obj<Mesh>& mesh, const Obj<Section>& partition, const int height = 0) {
1126: MeshManager<Mesh> manager(mesh, height);
1127: int *start = NULL;
1128: int *adjacency = NULL;
1130: if (height == 0) {
1131: int numVertices;
1133: buildDualCSRV(mesh, &numVertices, &start, &adjacency, GraphPartitioner::zeroBase());
1134: GraphPartitioner().partition(numVertices, start, adjacency, partition, manager);
1135: destroyCSR(numVertices, start, adjacency);
1136: } else if (height == 1) {
1137: int numEdges;
1139: throw ALE::Exception("Not yet implemented");
1140: #if 0
1141: buildFaceDualCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
1142: #endif
1143: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
1144: destroyCSR(numEdges, start, adjacency);
1145: } else {
1146: throw ALE::Exception("Invalid partition height");
1147: }
1148: };
1149: // Add in the points in the closure (and star) of the partitioned points
1150: template<typename Mesh, typename Section>
1151: static void createPartitionClosure(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1152: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1153: const typename Section::chart_type& chart = pointPartition->getChart();
1154: size_t size = 0;
1156: for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1157: const typename Section::value_type *points = pointPartition->restrictPoint(*r_iter);
1158: const int numPoints = pointPartition->getFiberDimension(*r_iter);
1159: std::set<typename Section::value_type> closure;
1161: // TODO: Use Quiver's closure() here instead
1162: for(int p = 0; p < numPoints; ++p) {
1163: Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1164: Obj<typename Mesh::sieve_type::coneSet> next = new typename Mesh::sieve_type::coneSet();
1165: Obj<typename Mesh::sieve_type::coneSet> tmp;
1167: current->insert(points[p]);
1168: closure.insert(points[p]);
1169: while(current->size()) {
1170: for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1171: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1172:
1173: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1174: closure.insert(*c_iter);
1175: next->insert(*c_iter);
1176: }
1177: }
1178: tmp = current; current = next; next = tmp;
1179: next->clear();
1180: }
1181: if (height) {
1182: current->insert(points[p]);
1183: while(current->size()) {
1184: for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1185: const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1186:
1187: for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1188: closure.insert(*s_iter);
1189: next->insert(*s_iter);
1190: }
1191: }
1192: tmp = current; current = next; next = tmp;
1193: next->clear();
1194: }
1195: }
1196: }
1197: partition->setFiberDimension(*r_iter, closure.size());
1198: size = std::max(size, closure.size());
1199: }
1200: partition->allocatePoint();
1201: typename Section::value_type *values = new typename Section::value_type[size];
1203: for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1204: const typename Section::value_type *points = pointPartition->restrictPoint(*r_iter);
1205: const int numPoints = pointPartition->getFiberDimension(*r_iter);
1206: std::set<typename Section::value_type> closure;
1208: // TODO: Use Quiver's closure() here instead
1209: for(int p = 0; p < numPoints; ++p) {
1210: Obj<typename Mesh::sieve_type::coneSet> current = new typename Mesh::sieve_type::coneSet();
1211: Obj<typename Mesh::sieve_type::coneSet> next = new typename Mesh::sieve_type::coneSet();
1212: Obj<typename Mesh::sieve_type::coneSet> tmp;
1214: current->insert(points[p]);
1215: closure.insert(points[p]);
1216: while(current->size()) {
1217: for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1218: const Obj<typename Mesh::sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
1219:
1220: for(typename Mesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
1221: closure.insert(*c_iter);
1222: next->insert(*c_iter);
1223: }
1224: }
1225: tmp = current; current = next; next = tmp;
1226: next->clear();
1227: }
1228: if (height) {
1229: current->insert(points[p]);
1230: while(current->size()) {
1231: for(typename Mesh::sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
1232: const Obj<typename Mesh::sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
1233:
1234: for(typename Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
1235: closure.insert(*s_iter);
1236: next->insert(*s_iter);
1237: }
1238: }
1239: tmp = current; current = next; next = tmp;
1240: next->clear();
1241: }
1242: }
1243: }
1244: int i = 0;
1246: for(typename std::set<typename Section::value_type>::const_iterator p_iter = closure.begin(); p_iter != closure.end(); ++p_iter, ++i) {
1247: values[i] = *p_iter;
1248: }
1249: partition->updatePoint(*r_iter, values);
1250: }
1251: delete [] values;
1252: };
1253: template<typename Mesh, typename Section>
1254: static void createPartitionClosureV(const Obj<Mesh>& mesh, const Obj<Section>& pointPartition, const Obj<Section>& partition, const int height = 0) {
1255: typedef ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type> visitor_type;
1256: const Obj<typename Mesh::sieve_type>& sieve = mesh->getSieve();
1257: const typename Section::chart_type& chart = pointPartition->getChart();
1258: size_t size = 0;
1260: for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1261: const typename Section::value_type *points = pointPartition->restrictPoint(*r_iter);
1262: const int numPoints = pointPartition->getFiberDimension(*r_iter);
1263: typename visitor_type::visitor_type nV;
1264: visitor_type cV(*sieve, nV);
1266: for(int p = 0; p < numPoints; ++p) {
1267: sieve->cone(points[p], cV);
1268: if (height) {
1269: cV.setIsCone(false);
1270: sieve->support(points[p], cV);
1271: }
1272: }
1273: partition->setFiberDimension(*r_iter, cV.getPoints().size());
1274: size = std::max(size, cV.getPoints().size());
1275: }
1276: partition->allocatePoint();
1277: typename Section::value_type *values = new typename Section::value_type[size];
1279: for(typename Section::chart_type::const_iterator r_iter = chart.begin(); r_iter != chart.end(); ++r_iter) {
1280: const typename Section::value_type *points = pointPartition->restrictPoint(*r_iter);
1281: const int numPoints = pointPartition->getFiberDimension(*r_iter);
1282: typename visitor_type::visitor_type nV;
1283: visitor_type cV(*sieve, nV);
1285: for(int p = 0; p < numPoints; ++p) {
1286: sieve->cone(points[p], cV);
1287: if (height) {
1288: cV.setIsCone(false);
1289: sieve->support(points[p], cV);
1290: }
1291: }
1292: int i = 0;
1294: for(typename std::set<typename Mesh::point_type>::const_iterator p_iter = cV.getPoints().begin(); p_iter != cV.getPoints().end(); ++p_iter, ++i) {
1295: values[i] = *p_iter;
1296: }
1297: partition->updatePoint(*r_iter, values);
1298: }
1299: delete [] values;
1300: };
1301: // Create a section mapping points to partitions
1302: template<typename Section, typename MapSection>
1303: static void createPartitionMap(const Obj<Section>& partition, const Obj<MapSection>& partitionMap) {
1304: const typename Section::chart_type& chart = partition->getChart();
1306: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1307: partitionMap->setFiberDimension(*p_iter, 1);
1308: }
1309: partitionMap->allocatePoint();
1310: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1311: const typename Section::value_type *points = partition->restrictPoint(*p_iter);
1312: const int size = partition->getFiberDimension(*p_iter);
1313: const typename Section::point_type part = *p_iter;
1315: for(int i = 0; i < size; ++i) {
1316: partitionMap->updatePoint(points[i], &part);
1317: }
1318: }
1319: };
1320: };
1321: #endif
1323: namespace New {
1324: template<typename Bundle_, typename Alloc_ = typename Bundle_::alloc_type>
1325: class Partitioner {
1326: public:
1327: typedef Bundle_ bundle_type;
1328: typedef Alloc_ alloc_type;
1329: typedef typename bundle_type::sieve_type sieve_type;
1330: typedef typename bundle_type::point_type point_type;
1331: public:
1334: // This creates a CSR representation of the adjacency matrix for cells
1335: // - We allow an exception to contiguous numbering.
1336: // If the cell id > numElements, we assign a new number starting at
1337: // the top and going downward. I know these might not match up with
1338: // the iterator order, but we can fix it later.
1339: static void buildDualCSR(const Obj<bundle_type>& bundle, const int dim, int **offsets, int **adjacency) {
1340: ALE_LOG_EVENT_BEGIN;
1341: typedef typename ALE::New::Completion<bundle_type, point_type, alloc_type> completion;
1342: const Obj<sieve_type>& sieve = bundle->getSieve();
1343: const Obj<typename bundle_type::label_sequence>& elements = bundle->heightStratum(0);
1344: Obj<sieve_type> overlapSieve = new sieve_type(bundle->comm(), bundle->debug());
1345: std::map<point_type, point_type> newCells;
1346: int numElements = elements->size();
1347: int newCell = numElements;
1348: int *off = new int[numElements+1];
1349: int offset = 0;
1350: int *adj;
1352: completion::scatterSupports(sieve, overlapSieve, bundle->getSendOverlap(), bundle->getRecvOverlap(), bundle);
1353: if (numElements == 0) {
1354: *offsets = NULL;
1355: *adjacency = NULL;
1356: ALE_LOG_EVENT_END;
1357: return;
1358: }
1359: if (bundle->depth() == dim) {
1360: int e = 1;
1362: off[0] = 0;
1363: for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1364: const Obj<typename sieve_type::traits::coneSequence>& faces = sieve->cone(*e_iter);
1365: typename sieve_type::traits::coneSequence::iterator fBegin = faces->begin();
1366: typename sieve_type::traits::coneSequence::iterator fEnd = faces->end();
1368: off[e] = off[e-1];
1369: for(typename sieve_type::traits::coneSequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1370: if (sieve->support(*f_iter)->size() == 2) {
1371: off[e]++;
1372: } else if ((sieve->support(*f_iter)->size() == 1) && (overlapSieve->support(*f_iter)->size() == 1)) {
1373: off[e]++;
1374: }
1375: }
1376: e++;
1377: }
1378: adj = new int[off[numElements]];
1379: for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1380: const Obj<typename sieve_type::traits::coneSequence>& faces = sieve->cone(*e_iter);
1381: typename sieve_type::traits::coneSequence::iterator fBegin = faces->begin();
1382: typename sieve_type::traits::coneSequence::iterator fEnd = faces->end();
1384: for(typename sieve_type::traits::coneSequence::iterator f_iter = fBegin; f_iter != fEnd; ++f_iter) {
1385: const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*f_iter);
1386: typename sieve_type::traits::supportSequence::iterator nBegin = neighbors->begin();
1387: typename sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
1389: for(typename sieve_type::traits::supportSequence::iterator n_iter = nBegin; n_iter != nEnd; ++n_iter) {
1390: if (*n_iter != *e_iter) adj[offset++] = *n_iter;
1391: }
1392: const Obj<typename sieve_type::traits::supportSequence>& oNeighbors = overlapSieve->support(*f_iter);
1393: typename sieve_type::traits::supportSequence::iterator onBegin = oNeighbors->begin();
1394: typename sieve_type::traits::supportSequence::iterator onEnd = oNeighbors->end();
1396: for(typename sieve_type::traits::supportSequence::iterator n_iter = onBegin; n_iter != onEnd; ++n_iter) {
1397: adj[offset++] = *n_iter;
1398: }
1399: }
1400: }
1401: } else if (bundle->depth() == 1) {
1402: std::set<point_type> *neighborCells = new std::set<point_type>[numElements];
1403: int corners = sieve->cone(*elements->begin())->size();
1404: int faceVertices = -1;
1406: if (corners == dim+1) {
1407: faceVertices = dim;
1408: } else if ((dim == 2) && (corners == 4)) {
1409: faceVertices = 2;
1410: } else if ((dim == 3) && (corners == 8)) {
1411: faceVertices = 4;
1412: } else {
1413: throw ALE::Exception("Could not determine number of face vertices");
1414: }
1415: for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1416: const Obj<typename sieve_type::traits::coneSequence>& vertices = sieve->cone(*e_iter);
1417: typename sieve_type::traits::coneSequence::iterator vEnd = vertices->end();
1419: for(typename sieve_type::traits::coneSequence::iterator v_iter = vertices->begin(); v_iter != vEnd; ++v_iter) {
1420: const Obj<typename sieve_type::traits::supportSequence>& neighbors = sieve->support(*v_iter);
1421: typename sieve_type::traits::supportSequence::iterator nEnd = neighbors->end();
1423: for(typename sieve_type::traits::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
1424: if (*e_iter == *n_iter) continue;
1425: if ((int) sieve->nMeet(*e_iter, *n_iter, 1)->size() == faceVertices) {
1426: if ((*e_iter < numElements) && (*n_iter < numElements)) {
1427: neighborCells[*e_iter].insert(*n_iter);
1428: } else {
1429: point_type e = *e_iter, n = *n_iter;
1431: if (*e_iter >= numElements) {
1432: if (newCells.find(*e_iter) == newCells.end()) newCells[*e_iter] = --newCell;
1433: e = newCells[*e_iter];
1434: }
1435: if (*n_iter >= numElements) {
1436: if (newCells.find(*n_iter) == newCells.end()) newCells[*n_iter] = --newCell;
1437: n = newCells[*n_iter];
1438: }
1439: neighborCells[e].insert(n);
1440: }
1441: }
1442: }
1443: }
1444: }
1445: off[0] = 0;
1446: for(int e = 1; e <= numElements; e++) {
1447: off[e] = neighborCells[e-1].size() + off[e-1];
1448: }
1449: adj = new int[off[numElements]];
1450: for(int e = 0; e < numElements; e++) {
1451: for(typename std::set<point_type>::iterator n_iter = neighborCells[e].begin(); n_iter != neighborCells[e].end(); ++n_iter) {
1452: adj[offset++] = *n_iter;
1453: }
1454: }
1455: delete [] neighborCells;
1456: } else {
1457: throw ALE::Exception("Dual creation not defined for partially interpolated meshes");
1458: }
1459: if (offset != off[numElements]) {
1460: ostringstream msg;
1461: msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1462: throw ALE::Exception(msg.str().c_str());
1463: }
1464: //std::cout << "numElements: " << numElements << " newCell: " << newCell << std::endl;
1465: *offsets = off;
1466: *adjacency = adj;
1467: ALE_LOG_EVENT_END;
1468: };
1471: // This creates a CSR representation of the adjacency hypergraph for faces
1472: 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) {
1473: ALE_LOG_EVENT_BEGIN;
1474: const Obj<sieve_type>& sieve = bundle->getSieve();
1475: const Obj<typename bundle_type::label_sequence>& elements = bundle->heightStratum(0);
1476: int numElements = elements->size();
1477: int *off = new int[numElements+1];
1478: int e;
1480: if (bundle->depth() != dim) {
1481: throw ALE::Exception("Not yet implemented for non-interpolated meshes");
1482: }
1483: off[0] = 0;
1484: e = 1;
1485: for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1486: off[e] = sieve->cone(*e_iter)->size() + off[e-1];
1487: e++;
1488: }
1489: int *adj = new int[off[numElements]];
1490: int offset = 0;
1491: for(typename bundle_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1492: const Obj<typename sieve_type::traits::coneSequence>& faces = sieve->cone(*e_iter);
1493: typename sieve_type::traits::coneSequence::iterator fEnd = faces->end();
1495: for(typename sieve_type::traits::coneSequence::iterator f_iter = faces->begin(); f_iter != fEnd; ++f_iter) {
1496: adj[offset++] = fNumbering->getIndex(*f_iter);
1497: }
1498: }
1499: if (offset != off[numElements]) {
1500: ostringstream msg;
1501: msg << "ERROR: Total number of neighbors " << offset << " does not match the offset array " << off[numElements];
1502: throw ALE::Exception(msg.str().c_str());
1503: }
1504: *numEdges = numElements;
1505: *offsets = off;
1506: *adjacency = adj;
1507: ALE_LOG_EVENT_END;
1508: };
1509: template<typename PartitionType>
1510: static PartitionType *subordinatePartition(const Obj<bundle_type>& bundle, int levels, const Obj<bundle_type>& subBundle, const PartitionType assignment[]) {
1511: const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth());
1512: const Obj<typename bundle_type::label_sequence>& cells = subBundle->heightStratum(0);
1513: const Obj<typename bundle_type::numbering_type>& sNumbering = bundle->getFactory()->getLocalNumbering(subBundle, subBundle->depth());
1514: const int numCells = cells->size();
1515: PartitionType *subAssignment = new PartitionType[numCells];
1517: if (levels != 1) {
1518: throw ALE::Exception("Cannot calculate subordinate partition for any level separation other than 1");
1519: } else {
1520: const Obj<typename bundle_type::sieve_type>& sieve = bundle->getSieve();
1521: const Obj<typename bundle_type::sieve_type>& subSieve = subBundle->getSieve();
1522: Obj<typename bundle_type::sieve_type::coneSet> tmpSet = new typename bundle_type::sieve_type::coneSet();
1523: Obj<typename bundle_type::sieve_type::coneSet> tmpSet2 = new typename bundle_type::sieve_type::coneSet();
1525: for(typename bundle_type::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
1526: const Obj<typename bundle_type::sieve_type::coneSequence>& cone = subSieve->cone(*c_iter);
1528: Obj<typename bundle_type::sieve_type::supportSet> cell = sieve->nJoin1(cone);
1529: if (cell->size() != 1) {
1530: std::cout << "Indeterminate subordinate partition for face " << *c_iter << std::endl;
1531: for(typename bundle_type::sieve_type::supportSet::iterator s_iter = cell->begin(); s_iter != cell->end(); ++s_iter) {
1532: std::cout << " cell " << *s_iter << std::endl;
1533: }
1534: // Could relax this to choosing the first one
1535: throw ALE::Exception("Indeterminate subordinate partition");
1536: }
1537: subAssignment[sNumbering->getIndex(*c_iter)] = assignment[cNumbering->getIndex(*cell->begin())];
1538: tmpSet->clear();
1539: tmpSet2->clear();
1540: }
1541: }
1542: return subAssignment;
1543: };
1544: };
1545: #ifdef PETSC_HAVE_CHACO
1546: namespace Chaco {
1547: template<typename Bundle_>
1548: class Partitioner {
1549: public:
1550: typedef Bundle_ bundle_type;
1551: typedef typename bundle_type::sieve_type sieve_type;
1552: typedef typename bundle_type::point_type point_type;
1553: typedef short int part_type;
1554: public:
1557: static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1558: part_type *assignment = NULL; /* set number of each vtx (length n) */
1559: int *start; /* start of edge list for each vertex */
1560: int *adjacency; /* = adj -> j; edge list data */
1562: ALE_LOG_EVENT_BEGIN;
1563: ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &start, &adjacency);
1564: if (bundle->commRank() == 0) {
1565: /* arguments for Chaco library */
1566: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
1567: int nvtxs; /* number of vertices in full graph */
1568: int *vwgts = NULL; /* weights for all vertices */
1569: float *ewgts = NULL; /* weights for all edges */
1570: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
1571: char *outassignname = NULL; /* name of assignment output file */
1572: char *outfilename = NULL; /* output file name */
1573: int architecture = dim; /* 0 => hypercube, d => d-dimensional mesh */
1574: int ndims_tot = 0; /* total number of cube dimensions to divide */
1575: int mesh_dims[3]; /* dimensions of mesh of processors */
1576: double *goal = NULL; /* desired set sizes for each set */
1577: int global_method = 1; /* global partitioning algorithm */
1578: int local_method = 1; /* local partitioning algorithm */
1579: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
1580: int vmax = 200; /* how many vertices to coarsen down to? */
1581: int ndims = 1; /* number of eigenvectors (2^d sets) */
1582: double eigtol = 0.001; /* tolerance on eigenvectors */
1583: long seed = 123636512; /* for random graph mutations */
1584: float *vCoords[3];
1587: PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_global_method", &global_method, PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1588: PetscOptionsGetInt(PETSC_NULL, "-partitioner_chaco_local_method", &local_method, PETSC_NULL);CHKERROR(ierr, "Error in PetscOptionsGetInt");
1589: if (global_method == 3) {
1590: // Inertial Partitioning
1591: PetscMalloc3(nvtxs,float,&x,nvtxs,float,&y,nvtxs,float,&z);CHKERROR(ierr, "Error in PetscMalloc");
1592: vCoords[0] = x; vCoords[1] = y; vCoords[2] = z;
1593: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(0);
1594: const Obj<typename bundle_type::real_section_type>& coordinates = bundle->getRealSection("coordinates");
1595: const int corners = bundle->size(coordinates, *(cells->begin()))/dim;
1596: int c = 0;
1598: for(typename bundle_type::label_sequence::iterator c_iter = cells->begin(); c_iter !=cells->end(); ++c_iter, ++c) {
1599: const double *coords = bundle->restrictClosure(coordinates, *c_iter);
1601: for(int d = 0; d < dim; ++d) {
1602: vCoords[d][c] = 0.0;
1603: }
1604: for(int v = 0; v < corners; ++v) {
1605: for(int d = 0; d < dim; ++d) {
1606: vCoords[d][c] += coords[v*dim+d];
1607: }
1608: }
1609: for(int d = 0; d < dim; ++d) {
1610: vCoords[d][c] /= corners;
1611: }
1612: }
1613: }
1615: nvtxs = bundle->heightStratum(0)->size();
1616: mesh_dims[0] = bundle->commSize(); mesh_dims[1] = 1; mesh_dims[2] = 1;
1617: for(int e = 0; e < start[nvtxs]; e++) {
1618: adjacency[e]++;
1619: }
1620: assignment = new part_type[nvtxs];
1621: PetscMemzero(assignment, nvtxs * sizeof(part_type));CHKERROR(ierr, "Error in PetscMemzero");
1623: /* redirect output to buffer: chaco -> msgLog */
1624: #ifdef PETSC_HAVE_UNISTD_H
1625: char *msgLog;
1626: int fd_stdout, fd_pipe[2], count;
1628: fd_stdout = dup(1);
1629: pipe(fd_pipe);
1630: close(1);
1631: dup2(fd_pipe[1], 1);
1632: msgLog = new char[16284];
1633: #endif
1635: interface(nvtxs, start, adjacency, vwgts, ewgts, x, y, z,
1636: outassignname, outfilename, assignment, architecture, ndims_tot,
1637: mesh_dims, goal, global_method, local_method, rqi_flag, vmax, ndims,
1638: eigtol, seed);
1640: #ifdef PETSC_HAVE_UNISTD_H
1641: int SIZE_LOG = 10000;
1643: fflush(stdout);
1644: count = read(fd_pipe[0], msgLog, (SIZE_LOG - 1) * sizeof(char));
1645: if (count < 0) count = 0;
1646: msgLog[count] = 0;
1647: close(1);
1648: dup2(fd_stdout, 1);
1649: close(fd_stdout);
1650: close(fd_pipe[0]);
1651: close(fd_pipe[1]);
1652: if (bundle->debug()) {
1653: std::cout << msgLog << std::endl;
1654: }
1655: delete [] msgLog;
1656: #endif
1657: if (global_method == 3) {
1658: // Inertial Partitioning
1659: PetscFree3(x, y, z);CHKERROR(ierr, "Error in PetscFree");
1660: }
1661: }
1662: if (adjacency) delete [] adjacency;
1663: if (start) delete [] start;
1664: ALE_LOG_EVENT_END;
1665: return assignment;
1666: };
1667: static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1668: throw ALE::Exception("Chaco cannot partition a mesh by faces");
1669: };
1670: };
1671: };
1672: #endif
1673: #ifdef PETSC_HAVE_PARMETIS
1674: namespace ParMetis {
1675: template<typename Bundle_>
1676: class Partitioner {
1677: public:
1678: typedef Bundle_ bundle_type;
1679: typedef typename bundle_type::sieve_type sieve_type;
1680: typedef typename bundle_type::point_type point_type;
1681: typedef int part_type;
1682: public:
1685: static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1686: int nvtxs = 0; // The number of vertices in full graph
1687: int *vtxdist; // Distribution of vertices across processes
1688: int *xadj; // Start of edge list for each vertex
1689: int *adjncy; // Edge lists for all vertices
1690: int *vwgt = NULL; // Vertex weights
1691: int *adjwgt = NULL; // Edge weights
1692: int wgtflag = 0; // Indicates which weights are present
1693: int numflag = 0; // Indicates initial offset (0 or 1)
1694: int ncon = 1; // The number of weights per vertex
1695: int nparts = bundle->commSize(); // The number of partitions
1696: float *tpwgts; // The fraction of vertex weights assigned to each partition
1697: float *ubvec; // The balance intolerance for vertex weights
1698: int options[5]; // Options
1699: // Outputs
1700: int edgeCut; // The number of edges cut by the partition
1701: int *assignment = NULL; // The vertex partition
1703: options[0] = 0; // Use all defaults
1704: vtxdist = new int[nparts+1];
1705: vtxdist[0] = 0;
1706: tpwgts = new float[ncon*nparts];
1707: for(int p = 0; p < nparts; ++p) {
1708: tpwgts[p] = 1.0/nparts;
1709: }
1710: ubvec = new float[ncon];
1711: ubvec[0] = 1.05;
1712: nvtxs = bundle->heightStratum(0)->size();
1713: assignment = new part_type[nvtxs];
1714: MPI_Allgather(&nvtxs, 1, MPI_INT, &vtxdist[1], 1, MPI_INT, bundle->comm());
1715: for(int p = 2; p <= nparts; ++p) {
1716: vtxdist[p] += vtxdist[p-1];
1717: }
1718: if (bundle->commSize() == 1) {
1719: PetscMemzero(assignment, nvtxs * sizeof(part_type));
1720: } else {
1721: ALE::New::Partitioner<bundle_type>::buildDualCSR(bundle, dim, &xadj, &adjncy);
1723: if (bundle->debug() && nvtxs) {
1724: for(int p = 0; p <= nvtxs; ++p) {
1725: std::cout << "["<<bundle->commRank()<<"]xadj["<<p<<"] = " << xadj[p] << std::endl;
1726: }
1727: for(int i = 0; i < xadj[nvtxs]; ++i) {
1728: std::cout << "["<<bundle->commRank()<<"]adjncy["<<i<<"] = " << adjncy[i] << std::endl;
1729: }
1730: }
1731: if (vtxdist[1] == vtxdist[nparts]) {
1732: if (bundle->commRank() == 0) {
1733: METIS_PartGraphKway(&nvtxs, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &nparts, options, &edgeCut, assignment);
1734: if (bundle->debug()) {std::cout << "Metis: edgecut is " << edgeCut << std::endl;}
1735: }
1736: } else {
1737: MPI_Comm comm = bundle->comm();
1739: ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
1740: if (bundle->debug()) {std::cout << "ParMetis: edgecut is " << edgeCut << std::endl;}
1741: }
1742: if (xadj != NULL) delete [] xadj;
1743: if (adjncy != NULL) delete [] adjncy;
1744: }
1745: delete [] vtxdist;
1746: delete [] tpwgts;
1747: delete [] ubvec;
1748: return assignment;
1749: };
1752: static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1753: #ifdef PETSC_HAVE_HMETIS
1754: int *assignment = NULL; // The vertex partition
1755: int nvtxs; // The number of vertices
1756: int nhedges; // The number of hyperedges
1757: int *vwgts; // The vertex weights
1758: int *eptr; // The offsets of each hyperedge
1759: int *eind; // The vertices in each hyperedge, indexed by eptr
1760: int *hewgts; // The hyperedge weights
1761: int nparts; // The number of partitions
1762: int ubfactor; // The allowed load imbalance (1-50)
1763: int options[9]; // Options
1764: // Outputs
1765: int edgeCut; // The number of edges cut by the partition
1766: const Obj<ALE::Mesh::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);
1768: if (topology->commRank() == 0) {
1769: nvtxs = bundle->heightStratum(1)->size();
1770: vwgts = NULL;
1771: hewgts = NULL;
1772: nparts = bundle->commSize();
1773: ubfactor = 5;
1774: options[0] = 1; // Use all defaults
1775: options[1] = 10; // Number of bisections tested
1776: options[2] = 1; // Vertex grouping scheme
1777: options[3] = 1; // Objective function
1778: options[4] = 1; // V-cycle refinement
1779: options[5] = 0;
1780: options[6] = 0;
1781: options[7] = 1; // Random seed
1782: options[8] = 24; // Debugging level
1783: assignment = new part_type[nvtxs];
1785: if (bundle->commSize() == 1) {
1786: PetscMemzero(assignment, nvtxs * sizeof(part_type));
1787: } else {
1788: ALE::New::Partitioner<bundle_type>::buildFaceCSR(bundle, dim, fNumbering, &nhedges, &eptr, &eind);
1789: HMETIS_PartKway(nvtxs, nhedges, vwgts, eptr, eind, hewgts, nparts, ubfactor, options, assignment, &edgeCut);
1791: delete [] eptr;
1792: delete [] eind;
1793: }
1794: if (bundle->debug()) {for (int i = 0; i<nvtxs; i++) printf("[%d] %d\n", PetscGlobalRank, assignment[i]);}
1795: } else {
1796: assignment = NULL;
1797: }
1798: return assignment;
1799: #else
1800: throw ALE::Exception("hmetis partitioner is not available.");
1801: #endif
1802: };
1803: };
1804: };
1805: #endif
1806: #ifdef PETSC_HAVE_ZOLTAN
1807: namespace Zoltan {
1808: template<typename Bundle_>
1809: class Partitioner {
1810: public:
1811: typedef Bundle_ bundle_type;
1812: typedef typename bundle_type::sieve_type sieve_type;
1813: typedef typename bundle_type::point_type point_type;
1814: typedef int part_type;
1815: public:
1816: static part_type *partitionSieve(const Obj<bundle_type>& bundle, const int dim) {
1817: throw ALE::Exception("Zoltan partition by cells not implemented");
1818: };
1821: static part_type *partitionSieveByFace(const Obj<bundle_type>& bundle, const int dim) {
1822: // Outputs
1823: float version; // The library version
1824: int changed; // Did the partition change?
1825: int numGidEntries; // Number of array entries for a single global ID (1)
1826: int numLidEntries; // Number of array entries for a single local ID (1)
1827: int numImport; // The number of imported points
1828: ZOLTAN_ID_PTR import_global_ids; // The imported points
1829: ZOLTAN_ID_PTR import_local_ids; // The imported points
1830: int *import_procs; // The proc each point was imported from
1831: int *import_to_part; // The partition of each imported point
1832: int numExport; // The number of exported points
1833: ZOLTAN_ID_PTR export_global_ids; // The exported points
1834: ZOLTAN_ID_PTR export_local_ids; // The exported points
1835: int *export_procs; // The proc each point was exported to
1836: int *export_to_part; // The partition assignment of all local points
1837: int *assignment; // The partition assignment of all local points
1838: const Obj<typename bundle_type::numbering_type>& fNumbering = bundle->getFactory()->getNumbering(bundle, bundle->depth()-1);
1840: if (bundle->commSize() == 1) {
1841: PetscMemzero(assignment, bundle->heightStratum(1)->size() * sizeof(part_type));
1842: } else {
1843: if (bundle->commRank() == 0) {
1844: nvtxs_Zoltan = bundle->heightStratum(1)->size();
1845: ALE::New::Partitioner<bundle_type>::buildFaceCSR(bundle, dim, fNumbering, &nhedges_Zoltan, &eptr_Zoltan, &eind_Zoltan);
1846: assignment = new int[nvtxs_Zoltan];
1847: } else {
1848: nvtxs_Zoltan = bundle->heightStratum(1)->size();
1849: nhedges_Zoltan = 0;
1850: eptr_Zoltan = new int[1];
1851: eind_Zoltan = new int[1];
1852: eptr_Zoltan[0] = 0;
1853: assignment = NULL;
1854: }
1856: int Zoltan_Initialize(0, NULL, &version);
1857: struct Zoltan_Struct *zz = Zoltan_Create(bundle->comm());
1858: // General parameters
1859: Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
1860: Zoltan_Set_Param(zz, "LB_METHOD", "PHG");
1861: Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTITION");
1862: // PHG parameters
1863: Zoltan_Set_Param(zz, "PHG_OUTPUT_LEVEL", "2");
1864: Zoltan_Set_Param(zz, "PHG_EDGE_SIZE_THRESHOLD", "1.0"); // Do not throw out dense edges
1865: // Call backs
1866: Zoltan_Set_Num_Obj_Fn(zz, getNumVertices_Zoltan, NULL);
1867: Zoltan_Set_Obj_List_Fn(zz, getLocalElements_Zoltan, NULL);
1868: Zoltan_Set_HG_Size_CS_Fn(zz, getHgSizes_Zoltan, NULL);
1869: Zoltan_Set_HG_CS_Fn(zz, getHg_Zoltan, NULL);
1870: // Debugging
1871: //Zoltan_Generate_Files(zz, "zoltan.debug", 1, 0, 0, 1); // if using hypergraph callbacks
1873: Zoltan_LB_Partition(zz, &changed, &numGidEntries, &numLidEntries,
1874: &numImport, &import_global_ids, &import_local_ids, &import_procs, &import_to_part,
1875: &numExport, &export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1876: for(int v = 0; v < nvtxs_Zoltan; ++v) {
1877: assignment[v] = export_to_part[v];
1878: }
1879: Zoltan_LB_Free_Part(&import_global_ids, &import_local_ids, &import_procs, &import_to_part);
1880: Zoltan_LB_Free_Part(&export_global_ids, &export_local_ids, &export_procs, &export_to_part);
1881: Zoltan_Destroy(&zz);
1883: delete [] eptr_Zoltan;
1884: delete [] eind_Zoltan;
1885: }
1886: if (assignment) {for (int i=0; i<nvtxs_Zoltan; i++) printf("[%d] %d\n",PetscGlobalRank,assignment[i]);}
1887: return assignment;
1888: };
1889: };
1890: };
1891: #endif
1892: }
1893: }
1894: #endif