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