Actual source code: plexpartition.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
6: {
7: const PetscInt maxFaceCases = 30;
8: PetscInt numFaceCases = 0;
9: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
10: PetscInt *off, *adj;
11: PetscInt *neighborCells = NULL;
12: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
16: /* For parallel partitioning, I think you have to communicate supports */
17: DMPlexGetDimension(dm, &dim);
18: cellDim = dim - cellHeight;
19: DMPlexGetDepth(dm, &depth);
20: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
21: if (cEnd - cStart == 0) {
22: if (numVertices) *numVertices = 0;
23: if (offsets) *offsets = NULL;
24: if (adjacency) *adjacency = NULL;
25: return(0);
26: }
27: numCells = cEnd - cStart;
28: faceDepth = depth - cellHeight;
29: if (dim == depth) {
30: PetscInt f, fStart, fEnd;
32: PetscCalloc1(numCells+1, &off);
33: /* Count neighboring cells */
34: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
35: for (f = fStart; f < fEnd; ++f) {
36: const PetscInt *support;
37: PetscInt supportSize;
38: DMPlexGetSupportSize(dm, f, &supportSize);
39: DMPlexGetSupport(dm, f, &support);
40: if (supportSize == 2) {
41: ++off[support[0]-cStart+1];
42: ++off[support[1]-cStart+1];
43: }
44: }
45: /* Prefix sum */
46: for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
47: if (adjacency) {
48: PetscInt *tmp;
50: PetscMalloc1(off[numCells], &adj);
51: PetscMalloc1((numCells+1), &tmp);
52: PetscMemcpy(tmp, off, (numCells+1) * sizeof(PetscInt));
53: /* Get neighboring cells */
54: for (f = fStart; f < fEnd; ++f) {
55: const PetscInt *support;
56: PetscInt supportSize;
57: DMPlexGetSupportSize(dm, f, &supportSize);
58: DMPlexGetSupport(dm, f, &support);
59: if (supportSize == 2) {
60: adj[tmp[support[0]-cStart]++] = support[1];
61: adj[tmp[support[1]-cStart]++] = support[0];
62: }
63: }
64: #if defined(PETSC_USE_DEBUG)
65: for (c = 0; c < cEnd-cStart; ++c) if (tmp[c] != off[c+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
66: #endif
67: PetscFree(tmp);
68: }
69: if (numVertices) *numVertices = numCells;
70: if (offsets) *offsets = off;
71: if (adjacency) *adjacency = adj;
72: return(0);
73: }
74: /* Setup face recognition */
75: if (faceDepth == 1) {
76: PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */
78: for (c = cStart; c < cEnd; ++c) {
79: PetscInt corners;
81: DMPlexGetConeSize(dm, c, &corners);
82: if (!cornersSeen[corners]) {
83: PetscInt nFV;
85: if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
86: cornersSeen[corners] = 1;
88: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
90: numFaceVertices[numFaceCases++] = nFV;
91: }
92: }
93: }
94: PetscCalloc1(numCells+1, &off);
95: /* Count neighboring cells */
96: for (cell = cStart; cell < cEnd; ++cell) {
97: PetscInt numNeighbors = PETSC_DETERMINE, n;
99: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, &numNeighbors, &neighborCells);
100: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
101: for (n = 0; n < numNeighbors; ++n) {
102: PetscInt cellPair[2];
103: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
104: PetscInt meetSize = 0;
105: const PetscInt *meet = NULL;
107: cellPair[0] = cell; cellPair[1] = neighborCells[n];
108: if (cellPair[0] == cellPair[1]) continue;
109: if (!found) {
110: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
111: if (meetSize) {
112: PetscInt f;
114: for (f = 0; f < numFaceCases; ++f) {
115: if (numFaceVertices[f] == meetSize) {
116: found = PETSC_TRUE;
117: break;
118: }
119: }
120: }
121: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
122: }
123: if (found) ++off[cell-cStart+1];
124: }
125: }
126: /* Prefix sum */
127: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
129: if (adjacency) {
130: PetscMalloc1(off[numCells], &adj);
131: /* Get neighboring cells */
132: for (cell = cStart; cell < cEnd; ++cell) {
133: PetscInt numNeighbors = PETSC_DETERMINE, n;
134: PetscInt cellOffset = 0;
136: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, &numNeighbors, &neighborCells);
137: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
138: for (n = 0; n < numNeighbors; ++n) {
139: PetscInt cellPair[2];
140: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
141: PetscInt meetSize = 0;
142: const PetscInt *meet = NULL;
144: cellPair[0] = cell; cellPair[1] = neighborCells[n];
145: if (cellPair[0] == cellPair[1]) continue;
146: if (!found) {
147: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
148: if (meetSize) {
149: PetscInt f;
151: for (f = 0; f < numFaceCases; ++f) {
152: if (numFaceVertices[f] == meetSize) {
153: found = PETSC_TRUE;
154: break;
155: }
156: }
157: }
158: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
159: }
160: if (found) {
161: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
162: ++cellOffset;
163: }
164: }
165: }
166: }
167: PetscFree(neighborCells);
168: if (numVertices) *numVertices = numCells;
169: if (offsets) *offsets = off;
170: if (adjacency) *adjacency = adj;
171: return(0);
172: }
174: #if defined(PETSC_HAVE_CHACO)
175: #if defined(PETSC_HAVE_UNISTD_H)
176: #include <unistd.h>
177: #endif
178: /* Chaco does not have an include file */
179: PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
180: float *ewgts, float *x, float *y, float *z, char *outassignname,
181: char *outfilename, short *assignment, int architecture, int ndims_tot,
182: int mesh_dims[3], double *goal, int global_method, int local_method,
183: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
185: extern int FREE_GRAPH;
189: PetscErrorCode DMPlexPartition_Chaco(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
190: {
191: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
192: MPI_Comm comm;
193: int nvtxs = numVertices; /* number of vertices in full graph */
194: int *vwgts = NULL; /* weights for all vertices */
195: float *ewgts = NULL; /* weights for all edges */
196: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
197: char *outassignname = NULL; /* name of assignment output file */
198: char *outfilename = NULL; /* output file name */
199: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
200: int ndims_tot = 0; /* total number of cube dimensions to divide */
201: int mesh_dims[3]; /* dimensions of mesh of processors */
202: double *goal = NULL; /* desired set sizes for each set */
203: int global_method = 1; /* global partitioning algorithm */
204: int local_method = 1; /* local partitioning algorithm */
205: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
206: int vmax = 200; /* how many vertices to coarsen down to? */
207: int ndims = 1; /* number of eigenvectors (2^d sets) */
208: double eigtol = 0.001; /* tolerance on eigenvectors */
209: long seed = 123636512; /* for random graph mutations */
210: short int *assignment; /* Output partition */
211: int fd_stdout, fd_pipe[2];
212: PetscInt *points;
213: PetscMPIInt commSize;
214: int i, v, p;
218: PetscObjectGetComm((PetscObject)dm,&comm);
219: MPI_Comm_size(comm, &commSize);
220: if (!numVertices) {
221: PetscSectionCreate(comm, partSection);
222: PetscSectionSetChart(*partSection, 0, commSize);
223: PetscSectionSetUp(*partSection);
224: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
225: return(0);
226: }
227: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
228: for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];
230: if (global_method == INERTIAL_METHOD) {
231: /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
232: SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
233: }
234: mesh_dims[0] = commSize;
235: mesh_dims[1] = 1;
236: mesh_dims[2] = 1;
237: PetscMalloc1(nvtxs, &assignment);
238: /* Chaco outputs to stdout. We redirect this to a buffer. */
239: /* TODO: check error codes for UNIX calls */
240: #if defined(PETSC_HAVE_UNISTD_H)
241: {
242: int piperet;
243: piperet = pipe(fd_pipe);
244: if (piperet) SETERRQ(comm,PETSC_ERR_SYS,"Could not create pipe");
245: fd_stdout = dup(1);
246: close(1);
247: dup2(fd_pipe[1], 1);
248: }
249: #endif
250: interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
251: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
252: vmax, ndims, eigtol, seed);
253: #if defined(PETSC_HAVE_UNISTD_H)
254: {
255: char msgLog[10000];
256: int count;
258: fflush(stdout);
259: count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
260: if (count < 0) count = 0;
261: msgLog[count] = 0;
262: close(1);
263: dup2(fd_stdout, 1);
264: close(fd_stdout);
265: close(fd_pipe[0]);
266: close(fd_pipe[1]);
267: if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
268: }
269: #endif
270: /* Convert to PetscSection+IS */
271: PetscSectionCreate(comm, partSection);
272: PetscSectionSetChart(*partSection, 0, commSize);
273: for (v = 0; v < nvtxs; ++v) {
274: PetscSectionAddDof(*partSection, assignment[v], 1);
275: }
276: PetscSectionSetUp(*partSection);
277: PetscMalloc1(nvtxs, &points);
278: for (p = 0, i = 0; p < commSize; ++p) {
279: for (v = 0; v < nvtxs; ++v) {
280: if (assignment[v] == p) points[i++] = v;
281: }
282: }
283: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
284: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
285: if (global_method == INERTIAL_METHOD) {
286: /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
287: }
288: PetscFree(assignment);
289: for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
290: return(0);
291: }
292: #endif
294: #if defined(PETSC_HAVE_PARMETIS)
295: #include <parmetis.h>
299: PetscErrorCode DMPlexPartition_ParMetis(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
300: {
301: MPI_Comm comm;
302: PetscInt nvtxs = numVertices; /* The number of vertices in full graph */
303: PetscInt *vtxdist; /* Distribution of vertices across processes */
304: PetscInt *xadj = start; /* Start of edge list for each vertex */
305: PetscInt *adjncy = adjacency; /* Edge lists for all vertices */
306: PetscInt *vwgt = NULL; /* Vertex weights */
307: PetscInt *adjwgt = NULL; /* Edge weights */
308: PetscInt wgtflag = 0; /* Indicates which weights are present */
309: PetscInt numflag = 0; /* Indicates initial offset (0 or 1) */
310: PetscInt ncon = 1; /* The number of weights per vertex */
311: PetscInt nparts; /* The number of partitions */
312: PetscReal *tpwgts; /* The fraction of vertex weights assigned to each partition */
313: PetscReal *ubvec; /* The balance intolerance for vertex weights */
314: PetscInt options[5]; /* Options */
315: /* Outputs */
316: PetscInt edgeCut; /* The number of edges cut by the partition */
317: PetscInt *assignment, *points;
318: PetscMPIInt commSize, rank, p, v, i;
322: PetscObjectGetComm((PetscObject) dm, &comm);
323: MPI_Comm_size(comm, &commSize);
324: MPI_Comm_rank(comm, &rank);
325: nparts = commSize;
326: options[0] = 0; /* Use all defaults */
327: /* Calculate vertex distribution */
328: PetscMalloc4(nparts+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment);
329: vtxdist[0] = 0;
330: MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
331: for (p = 2; p <= nparts; ++p) {
332: vtxdist[p] += vtxdist[p-1];
333: }
334: /* Calculate weights */
335: for (p = 0; p < nparts; ++p) {
336: tpwgts[p] = 1.0/nparts;
337: }
338: ubvec[0] = 1.05;
340: if (nparts == 1) {
341: PetscMemzero(assignment, nvtxs * sizeof(PetscInt));
342: } else {
343: if (vtxdist[1] == vtxdist[nparts]) {
344: if (!rank) {
345: PetscStackPush("METIS_PartGraphKway");
346: METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, NULL, &edgeCut, assignment);
347: PetscStackPop;
348: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
349: }
350: } else {
351: PetscStackPush("ParMETIS_V3_PartKway");
352: ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
353: PetscStackPop;
354: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
355: }
356: }
357: /* Convert to PetscSection+IS */
358: PetscSectionCreate(comm, partSection);
359: PetscSectionSetChart(*partSection, 0, commSize);
360: for (v = 0; v < nvtxs; ++v) {
361: PetscSectionAddDof(*partSection, assignment[v], 1);
362: }
363: PetscSectionSetUp(*partSection);
364: PetscMalloc1(nvtxs, &points);
365: for (p = 0, i = 0; p < commSize; ++p) {
366: for (v = 0; v < nvtxs; ++v) {
367: if (assignment[v] == p) points[i++] = v;
368: }
369: }
370: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
371: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
372: PetscFree4(vtxdist,tpwgts,ubvec,assignment);
373: return(0);
374: }
375: #endif
379: /* Expand the partition by BFS on the adjacency graph */
380: PetscErrorCode DMPlexEnlargePartition(DM dm, const PetscInt start[], const PetscInt adjacency[], PetscSection origPartSection, IS origPartition, PetscSection *partSection, IS *partition)
381: {
382: PetscHashI h;
383: const PetscInt *points;
384: PetscInt **tmpPoints, *newPoints, totPoints = 0;
385: PetscInt pStart, pEnd, part, q;
386: PetscBool useCone;
387: PetscErrorCode ierr;
390: PetscHashICreate(h);
391: PetscSectionCreate(PetscObjectComm((PetscObject)dm), partSection);
392: PetscSectionGetChart(origPartSection, &pStart, &pEnd);
393: PetscSectionSetChart(*partSection, pStart, pEnd);
394: ISGetIndices(origPartition, &points);
395: PetscMalloc1((pEnd - pStart), &tmpPoints);
396: DMPlexGetAdjacencyUseCone(dm, &useCone);
397: DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);
398: for (part = pStart; part < pEnd; ++part) {
399: PetscInt *adj = NULL;
400: PetscInt numPoints, nP, numNewPoints, off, p, n = 0;
402: PetscHashIClear(h);
403: PetscSectionGetDof(origPartSection, part, &numPoints);
404: PetscSectionGetOffset(origPartSection, part, &off);
405: /* Add all existing points to h */
406: for (p = 0; p < numPoints; ++p) {
407: const PetscInt point = points[off+p];
408: PetscHashIAdd(h, point, 1);
409: }
410: PetscHashISize(h, nP);
411: if (nP != numPoints) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid partition has %d points, but only %d were unique", numPoints, nP);
412: /* Add all points in next BFS level */
413: for (p = 0; p < numPoints; ++p) {
414: const PetscInt point = points[off+p];
415: PetscInt adjSize = PETSC_DETERMINE, a;
417: DMPlexGetAdjacency(dm, point, &adjSize, &adj);
418: for (a = 0; a < adjSize; ++a) PetscHashIAdd(h, adj[a], 1);
419: }
420: PetscHashISize(h, numNewPoints);
421: PetscSectionSetDof(*partSection, part, numNewPoints);
422: PetscMalloc1(numNewPoints, &tmpPoints[part]);
423: PetscHashIGetKeys(h, &n, tmpPoints[part]);
424: PetscFree(adj);
425: totPoints += numNewPoints;
426: }
427: DMPlexSetAdjacencyUseCone(dm, useCone);
428: ISRestoreIndices(origPartition, &points);
429: PetscHashIDestroy(h);
430: PetscSectionSetUp(*partSection);
431: PetscMalloc1(totPoints, &newPoints);
432: for (part = pStart, q = 0; part < pEnd; ++part) {
433: PetscInt numPoints, p;
435: PetscSectionGetDof(*partSection, part, &numPoints);
436: for (p = 0; p < numPoints; ++p, ++q) newPoints[q] = tmpPoints[part][p];
437: PetscFree(tmpPoints[part]);
438: }
439: PetscFree(tmpPoints);
440: ISCreateGeneral(PetscObjectComm((PetscObject)dm), totPoints, newPoints, PETSC_OWN_POINTER, partition);
441: return(0);
442: }
446: /*
447: DMPlexCreatePartition - Create a non-overlapping partition of the points at the given height
449: Collective on DM
451: Input Parameters:
452: + dm - The DM
453: . height - The height for points in the partition
454: - enlarge - Expand each partition with neighbors
456: Output Parameters:
457: + partSection - The PetscSection giving the division of points by partition
458: . partition - The list of points by partition
459: . origPartSection - If enlarge is true, the PetscSection giving the division of points before enlarging by partition, otherwise NULL
460: - origPartition - If enlarge is true, the list of points before enlarging by partition, otherwise NULL
462: Level: developer
464: .seealso DMPlexDistribute()
465: */
466: PetscErrorCode DMPlexCreatePartition(DM dm, const char name[], PetscInt height, PetscBool enlarge, PetscSection *partSection, IS *partition, PetscSection *origPartSection, IS *origPartition)
467: {
468: char partname[1024];
469: PetscBool isChaco = PETSC_FALSE, isMetis = PETSC_FALSE, flg;
470: PetscMPIInt size;
474: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
476: *origPartSection = NULL;
477: *origPartition = NULL;
478: if (size == 1) {
479: PetscInt *points;
480: PetscInt cStart, cEnd, c;
482: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
483: PetscSectionCreate(PetscObjectComm((PetscObject)dm), partSection);
484: PetscSectionSetChart(*partSection, 0, size);
485: PetscSectionSetDof(*partSection, 0, cEnd-cStart);
486: PetscSectionSetUp(*partSection);
487: PetscMalloc1((cEnd - cStart), &points);
488: for (c = cStart; c < cEnd; ++c) points[c] = c;
489: ISCreateGeneral(PetscObjectComm((PetscObject)dm), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
490: return(0);
491: }
492: PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_partitioner", partname, 1024, &flg);
493: if (flg) name = partname;
494: if (name) {
495: PetscStrcmp(name, "chaco", &isChaco);
496: PetscStrcmp(name, "metis", &isMetis);
497: }
498: if (height == 0) {
499: PetscInt numVertices;
500: PetscInt *start = NULL;
501: PetscInt *adjacency = NULL;
503: DMPlexCreateNeighborCSR(dm, 0, &numVertices, &start, &adjacency);
504: if (!name || isChaco) {
505: #if defined(PETSC_HAVE_CHACO)
506: DMPlexPartition_Chaco(dm, numVertices, start, adjacency, partSection, partition);
507: #else
508: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
509: #endif
510: } else if (isMetis) {
511: #if defined(PETSC_HAVE_PARMETIS)
512: DMPlexPartition_ParMetis(dm, numVertices, start, adjacency, partSection, partition);
513: #endif
514: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unknown mesh partitioning package %s", name);
515: if (enlarge) {
516: *origPartSection = *partSection;
517: *origPartition = *partition;
519: DMPlexEnlargePartition(dm, start, adjacency, *origPartSection, *origPartition, partSection, partition);
520: }
521: PetscFree(start);
522: PetscFree(adjacency);
523: # if 0
524: } else if (height == 1) {
525: /* Build the dual graph for faces and partition the hypergraph */
526: PetscInt numEdges;
528: buildFaceCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
529: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
530: destroyCSR(numEdges, start, adjacency);
531: #endif
532: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid partition height %D", height);
533: return(0);
534: }
538: PetscErrorCode DMPlexCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition)
539: {
540: /* const PetscInt height = 0; */
541: const PetscInt *partArray;
542: PetscInt *allPoints, *packPoints;
543: PetscInt rStart, rEnd, rank, pStart, pEnd, newSize;
544: PetscErrorCode ierr;
545: PetscBT bt;
546: PetscSegBuffer segpack,segpart;
549: PetscSectionGetChart(pointSection, &rStart, &rEnd);
550: ISGetIndices(pointPartition, &partArray);
551: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
552: PetscSectionSetChart(*section, rStart, rEnd);
553: DMPlexGetChart(dm,&pStart,&pEnd);
554: PetscBTCreate(pEnd-pStart,&bt);
555: PetscSegBufferCreate(sizeof(PetscInt),1000,&segpack);
556: PetscSegBufferCreate(sizeof(PetscInt),1000,&segpart);
557: for (rank = rStart; rank < rEnd; ++rank) {
558: PetscInt partSize = 0, numPoints, offset, p, *PETSC_RESTRICT placePoints;
560: PetscSectionGetDof(pointSection, rank, &numPoints);
561: PetscSectionGetOffset(pointSection, rank, &offset);
562: for (p = 0; p < numPoints; ++p) {
563: PetscInt point = partArray[offset+p], closureSize, c;
564: PetscInt *closure = NULL;
566: /* TODO Include support for height > 0 case */
567: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
568: for (c=0; c<closureSize; c++) {
569: PetscInt cpoint = closure[c*2];
570: if (!PetscBTLookupSet(bt,cpoint-pStart)) {
571: PetscInt *PETSC_RESTRICT pt;
572: partSize++;
573: PetscSegBufferGetInts(segpart,1,&pt);
574: *pt = cpoint;
575: }
576: }
577: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
578: }
579: PetscSectionSetDof(*section, rank, partSize);
580: PetscSegBufferGetInts(segpack,partSize,&placePoints);
581: PetscSegBufferExtractTo(segpart,placePoints);
582: PetscSortInt(partSize,placePoints);
583: for (p=0; p<partSize; p++) {PetscBTClear(bt,placePoints[p]-pStart);}
584: }
585: PetscBTDestroy(&bt);
586: PetscSegBufferDestroy(&segpart);
588: PetscSectionSetUp(*section);
589: PetscSectionGetStorageSize(*section, &newSize);
590: PetscMalloc1(newSize, &allPoints);
592: PetscSegBufferExtractInPlace(segpack,&packPoints);
593: for (rank = rStart; rank < rEnd; ++rank) {
594: PetscInt numPoints, offset;
596: PetscSectionGetDof(*section, rank, &numPoints);
597: PetscSectionGetOffset(*section, rank, &offset);
598: PetscMemcpy(&allPoints[offset], packPoints, numPoints * sizeof(PetscInt));
599: packPoints += numPoints;
600: }
602: PetscSegBufferDestroy(&segpack);
603: ISRestoreIndices(pointPartition, &partArray);
604: ISCreateGeneral(PetscObjectComm((PetscObject)dm), newSize, allPoints, PETSC_OWN_POINTER, partition);
605: return(0);
606: }