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