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