1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/matorderimpl.h> /*I "petscmat.h" I*/
3 #include <petsc/private/dmlabelimpl.h>
4
DMPlexCreateOrderingClosure_Static(DM dm,PetscInt numPoints,const PetscInt pperm[],PetscInt ** clperm,PetscInt ** invclperm)5 static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
6 {
7 PetscInt *perm, *iperm;
8 PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p;
9
10 PetscFunctionBegin;
11 PetscCall(DMPlexGetDepth(dm, &depth));
12 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
13 PetscCall(PetscMalloc1(pEnd - pStart, &perm));
14 PetscCall(PetscMalloc1(pEnd - pStart, &iperm));
15 for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
16 for (d = depth; d > 0; --d) {
17 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
18 PetscCall(DMPlexGetDepthStratum(dm, d - 1, &fStart, &fEnd));
19 fMax = fStart;
20 for (p = pStart; p < pEnd; ++p) {
21 const PetscInt *cone;
22 PetscInt point, coneSize, c;
23
24 if (d == depth) {
25 perm[p] = pperm[p];
26 iperm[pperm[p]] = p;
27 }
28 point = perm[p];
29 PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
30 PetscCall(DMPlexGetCone(dm, point, &cone));
31 for (c = 0; c < coneSize; ++c) {
32 const PetscInt oldc = cone[c];
33 const PetscInt newc = iperm[oldc];
34
35 if (newc < 0) {
36 perm[fMax] = oldc;
37 iperm[oldc] = fMax++;
38 }
39 }
40 }
41 PetscCheck(fMax == fEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %" PetscInt_FMT " faces %" PetscInt_FMT " does not match permuted number %" PetscInt_FMT, d, fEnd - fStart, fMax - fStart);
42 }
43 *clperm = perm;
44 *invclperm = iperm;
45 PetscFunctionReturn(PETSC_SUCCESS);
46 }
47
48 /*@C
49 DMPlexGetOrdering - Calculate a reordering of the mesh
50
51 Collective
52
53 Input Parameters:
54 + dm - The `DMPLEX` object
55 . otype - type of reordering, see `MatOrderingType`
56 - label - [Optional] Label used to segregate ordering into sets, or `NULL`
57
58 Output Parameter:
59 . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
60
61 Level: intermediate
62
63 Note:
64 The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which
65 has different types of cells, and then loop over each set of reordered cells for assembly.
66
67 .seealso: `DMPLEX`, `DMPlexPermute()`, `MatOrderingType`, `MatGetOrdering()`
68 @*/
DMPlexGetOrdering(DM dm,MatOrderingType otype,DMLabel label,IS * perm)69 PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm)
70 {
71 PetscInt numCells = 0;
72 PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i;
73
74 PetscFunctionBegin;
75 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
76 PetscAssertPointer(perm, 4);
77 PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency));
78 PetscCall(PetscMalloc3(numCells, &cperm, numCells, &mask, numCells * 2, &xls));
79 if (numCells) {
80 /* Shift for Fortran numbering */
81 for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
82 for (i = 0; i <= numCells; ++i) ++start[i];
83 PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls));
84 }
85 PetscCall(PetscFree(start));
86 PetscCall(PetscFree(adjacency));
87 /* Shift for Fortran numbering */
88 for (c = 0; c < numCells; ++c) --cperm[c];
89 /* Segregate */
90 if (label) {
91 IS valueIS;
92 const PetscInt *valuesTmp;
93 PetscInt *values;
94 PetscInt numValues, numPoints = 0;
95 PetscInt *sperm, *vsize, *voff, v;
96
97 // Can't directly sort the valueIS, since it is a view into the DMLabel
98 PetscCall(DMLabelGetValueIS(label, &valueIS));
99 PetscCall(ISGetLocalSize(valueIS, &numValues));
100 PetscCall(ISGetIndices(valueIS, &valuesTmp));
101 PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff));
102 PetscCall(PetscArraycpy(values, valuesTmp, numValues));
103 PetscCall(PetscSortInt(numValues, values));
104 PetscCall(ISRestoreIndices(valueIS, &valuesTmp));
105 PetscCall(ISDestroy(&valueIS));
106 for (v = 0; v < numValues; ++v) {
107 PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v]));
108 if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1];
109 numPoints += vsize[v];
110 }
111 PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells);
112 for (c = 0; c < numCells; ++c) {
113 const PetscInt oldc = cperm[c];
114 PetscInt val, vloc;
115
116 PetscCall(DMLabelGetValue(label, oldc, &val));
117 PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc);
118 PetscCall(PetscFindInt(val, numValues, values, &vloc));
119 PetscCheck(vloc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val);
120 sperm[voff[vloc + 1]++] = oldc;
121 }
122 for (v = 0; v < numValues; ++v) PetscCheck(voff[v + 1] - voff[v] == vsize[v], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %" PetscInt_FMT " values found is %" PetscInt_FMT " != %" PetscInt_FMT, values[v], voff[v + 1] - voff[v], vsize[v]);
123 PetscCall(PetscArraycpy(cperm, sperm, numCells));
124 PetscCall(PetscFree4(sperm, values, vsize, voff));
125 }
126 /* Construct closure */
127 PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm));
128 PetscCall(PetscFree3(cperm, mask, xls));
129 PetscCall(PetscFree(clperm));
130 /* Invert permutation */
131 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
132 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
136 /*@
137 DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line
138
139 Collective
140
141 Input Parameter:
142 . dm - The `DMPLEX` object
143
144 Output Parameter:
145 . perm - The point permutation as an `IS`, `perm`[old point number] = new point number
146
147 Level: intermediate
148
149 .seealso: `DMPLEX`, `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()`
150 @*/
DMPlexGetOrdering1D(DM dm,IS * perm)151 PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm)
152 {
153 PetscInt *points;
154 const PetscInt *support, *cone;
155 PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0;
156
157 PetscFunctionBegin;
158 PetscCall(DMGetDimension(dm, &dim));
159 PetscCheck(dim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim);
160 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
161 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
162 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
163 PetscCall(PetscMalloc1(pEnd - pStart, &points));
164 for (c = cStart; c < cEnd; ++c) points[c] = c;
165 for (v = vStart; v < vEnd; ++v) points[v] = v;
166 for (v = vStart; v < vEnd; ++v) {
167 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
168 PetscCall(DMPlexGetSupport(dm, v, &support));
169 if (suppSize == 1) {
170 lastCell = support[0];
171 break;
172 }
173 }
174 if (v < vEnd) {
175 PetscInt pos = cEnd;
176
177 points[v] = pos++;
178 while (lastCell >= cStart) {
179 PetscCall(DMPlexGetCone(dm, lastCell, &cone));
180 if (cone[0] == v) v = cone[1];
181 else v = cone[0];
182 PetscCall(DMPlexGetSupport(dm, v, &support));
183 PetscCall(DMPlexGetSupportSize(dm, v, &suppSize));
184 if (suppSize == 1) {
185 lastCell = -1;
186 } else {
187 if (support[0] == lastCell) lastCell = support[1];
188 else lastCell = support[0];
189 }
190 points[v] = pos++;
191 }
192 PetscCheck(pos == pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd);
193 }
194 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm));
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
DMPlexRemapCoordinates_Private(IS perm,PetscSection cs,Vec coordinates,PetscSection * csNew,Vec * coordinatesNew)198 static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew)
199 {
200 PetscScalar *coords, *coordsNew;
201 const PetscInt *pperm;
202 PetscInt pStart, pEnd, p;
203 const char *name;
204
205 PetscFunctionBegin;
206 PetscCall(PetscSectionPermute(cs, perm, csNew));
207 PetscCall(VecDuplicate(coordinates, coordinatesNew));
208 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
209 PetscCall(PetscObjectSetName((PetscObject)*coordinatesNew, name));
210 PetscCall(VecGetArray(coordinates, &coords));
211 PetscCall(VecGetArray(*coordinatesNew, &coordsNew));
212 PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd));
213 PetscCall(ISGetIndices(perm, &pperm));
214 for (p = pStart; p < pEnd; ++p) {
215 PetscInt dof, off, offNew, d;
216
217 PetscCall(PetscSectionGetDof(*csNew, p, &dof));
218 PetscCall(PetscSectionGetOffset(cs, p, &off));
219 PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew));
220 for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d];
221 }
222 PetscCall(ISRestoreIndices(perm, &pperm));
223 PetscCall(VecRestoreArray(coordinates, &coords));
224 PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew));
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
228 /*@
229 DMPlexPermute - Reorder the mesh according to the input permutation
230
231 Collective
232
233 Input Parameters:
234 + dm - The `DMPLEX` object
235 - perm - The point permutation, `perm`[old point number] = new point number
236
237 Output Parameter:
238 . pdm - The permuted `DM`
239
240 Level: intermediate
241
242 .seealso: `DMPLEX`, `MatPermute()`
243 @*/
DMPlexPermute(DM dm,IS perm,DM * pdm)244 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
245 {
246 DM_Plex *plex = (DM_Plex *)dm->data, *plexNew;
247 PetscInt dim, cdim;
248 const char *name;
249
250 PetscFunctionBegin;
251 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
252 PetscValidHeaderSpecific(perm, IS_CLASSID, 2);
253 PetscAssertPointer(pdm, 3);
254 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), pdm));
255 PetscCall(DMSetType(*pdm, DMPLEX));
256 PetscCall(PetscObjectGetName((PetscObject)dm, &name));
257 PetscCall(PetscObjectSetName((PetscObject)*pdm, name));
258 PetscCall(DMGetDimension(dm, &dim));
259 PetscCall(DMSetDimension(*pdm, dim));
260 PetscCall(DMGetCoordinateDim(dm, &cdim));
261 PetscCall(DMSetCoordinateDim(*pdm, cdim));
262 PetscCall(DMCopyDisc(dm, *pdm));
263 if (dm->localSection) {
264 PetscSection section, sectionNew;
265
266 PetscCall(DMGetLocalSection(dm, §ion));
267 PetscCall(PetscSectionPermute(section, perm, §ionNew));
268 PetscCall(DMSetLocalSection(*pdm, sectionNew));
269 PetscCall(PetscSectionDestroy(§ionNew));
270 }
271 plexNew = (DM_Plex *)(*pdm)->data;
272 /* Ignore ltogmap, ltogmapb */
273 /* Ignore sf, sectionSF */
274 /* Ignore globalVertexNumbers, globalCellNumbers */
275 /* Reorder labels */
276 {
277 PetscInt numLabels, l;
278 DMLabel label, labelNew;
279
280 PetscCall(DMGetNumLabels(dm, &numLabels));
281 for (l = 0; l < numLabels; ++l) {
282 PetscCall(DMGetLabelByNum(dm, l, &label));
283 PetscCall(DMLabelPermute(label, perm, &labelNew));
284 PetscCall(DMAddLabel(*pdm, labelNew));
285 PetscCall(DMLabelDestroy(&labelNew));
286 }
287 PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel));
288 if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap));
289 }
290 if ((*pdm)->celltypeLabel) {
291 DMLabel ctLabel;
292
293 // Reset label for fast lookup
294 PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel));
295 PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel));
296 }
297 /* Reorder topology */
298 {
299 const PetscInt *pperm;
300 PetscInt n, pStart, pEnd, p;
301
302 PetscCall(PetscSectionDestroy(&plexNew->coneSection));
303 PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection));
304 PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n));
305 PetscCall(PetscMalloc1(n, &plexNew->cones));
306 PetscCall(PetscMalloc1(n, &plexNew->coneOrientations));
307 PetscCall(ISGetIndices(perm, &pperm));
308 PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd));
309 for (p = pStart; p < pEnd; ++p) {
310 PetscInt dof, off, offNew, d;
311
312 PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof));
313 PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off));
314 PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew));
315 for (d = 0; d < dof; ++d) {
316 plexNew->cones[offNew + d] = pperm[plex->cones[off + d]];
317 plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d];
318 }
319 }
320 PetscCall(PetscSectionDestroy(&plexNew->supportSection));
321 PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection));
322 PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n));
323 PetscCall(PetscMalloc1(n, &plexNew->supports));
324 PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd));
325 for (p = pStart; p < pEnd; ++p) {
326 PetscInt dof, off, offNew, d;
327
328 PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof));
329 PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off));
330 PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew));
331 for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]];
332 }
333 PetscCall(ISRestoreIndices(perm, &pperm));
334 }
335 /* Remap coordinates */
336 {
337 DM cdm, cdmNew;
338 PetscSection cs, csNew;
339 Vec coordinates, coordinatesNew;
340
341 PetscCall(DMGetCoordinateSection(dm, &cs));
342 PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
343 PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
344 PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
345 PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew));
346 PetscCall(PetscSectionDestroy(&csNew));
347 PetscCall(VecDestroy(&coordinatesNew));
348
349 PetscCall(DMGetCellCoordinateDM(dm, &cdm));
350 if (cdm) {
351 PetscCall(DMGetCoordinateDM(*pdm, &cdm));
352 PetscCall(DMClone(cdm, &cdmNew));
353 PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew));
354 PetscCall(DMDestroy(&cdmNew));
355 PetscCall(DMGetCellCoordinateSection(dm, &cs));
356 PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates));
357 PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew));
358 PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew));
359 PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew));
360 PetscCall(PetscSectionDestroy(&csNew));
361 PetscCall(VecDestroy(&coordinatesNew));
362 }
363 }
364 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm));
365 (*pdm)->setupcalled = PETSC_TRUE;
366 PetscFunctionReturn(PETSC_SUCCESS);
367 }
368
DMPlexReorderSetDefault_Plex(DM dm,DMReorderDefaultFlag reorder)369 PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
370 {
371 DM_Plex *mesh = (DM_Plex *)dm->data;
372
373 PetscFunctionBegin;
374 mesh->reorderDefault = reorder;
375 PetscFunctionReturn(PETSC_SUCCESS);
376 }
377
378 /*@
379 DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default
380
381 Logically Collective
382
383 Input Parameters:
384 + dm - The `DM`
385 - reorder - Flag for reordering
386
387 Level: intermediate
388
389 .seealso: `DMPlexReorderGetDefault()`
390 @*/
DMPlexReorderSetDefault(DM dm,DMReorderDefaultFlag reorder)391 PetscErrorCode DMPlexReorderSetDefault(DM dm, DMReorderDefaultFlag reorder)
392 {
393 PetscFunctionBegin;
394 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
395 PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMReorderDefaultFlag), (dm, reorder));
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
DMPlexReorderGetDefault_Plex(DM dm,DMReorderDefaultFlag * reorder)399 PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
400 {
401 DM_Plex *mesh = (DM_Plex *)dm->data;
402
403 PetscFunctionBegin;
404 *reorder = mesh->reorderDefault;
405 PetscFunctionReturn(PETSC_SUCCESS);
406 }
407
408 /*@
409 DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default
410
411 Not Collective
412
413 Input Parameter:
414 . dm - The `DM`
415
416 Output Parameter:
417 . reorder - Flag for reordering
418
419 Level: intermediate
420
421 .seealso: `DMPlexReorderSetDefault()`
422 @*/
DMPlexReorderGetDefault(DM dm,DMReorderDefaultFlag * reorder)423 PetscErrorCode DMPlexReorderGetDefault(DM dm, DMReorderDefaultFlag *reorder)
424 {
425 PetscFunctionBegin;
426 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
427 PetscAssertPointer(reorder, 2);
428 PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMReorderDefaultFlag *), (dm, reorder));
429 PetscFunctionReturn(PETSC_SUCCESS);
430 }
431
DMCreateSectionPermutation_Plex_Reverse(DM dm,IS * permutation,PetscBT * blockStarts)432 static PetscErrorCode DMCreateSectionPermutation_Plex_Reverse(DM dm, IS *permutation, PetscBT *blockStarts)
433 {
434 IS permIS;
435 PetscInt *perm;
436 PetscInt pStart, pEnd;
437
438 PetscFunctionBegin;
439 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
440 PetscCall(PetscMalloc1(pEnd - pStart, &perm));
441 for (PetscInt p = pStart; p < pEnd; ++p) perm[pEnd - 1 - p] = p;
442 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
443 PetscCall(ISSetPermutation(permIS));
444 *permutation = permIS;
445 PetscFunctionReturn(PETSC_SUCCESS);
446 }
447
448 // Reorder to group split nodes
DMCreateSectionPermutation_Plex_Cohesive_Old(DM dm,IS * permutation,PetscBT * blockStarts)449 static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive_Old(DM dm, IS *permutation, PetscBT *blockStarts)
450 {
451 IS permIS;
452 PetscBT bt, blst;
453 PetscInt *perm;
454 PetscInt pStart, pEnd, i = 0;
455
456 PetscFunctionBegin;
457 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
458 PetscCall(PetscMalloc1(pEnd - pStart, &perm));
459 PetscCall(PetscBTCreate(pEnd - pStart, &bt));
460 PetscCall(PetscBTCreate(pEnd - pStart, &blst));
461 for (PetscInt p = pStart; p < pEnd; ++p) {
462 const PetscInt *supp, *cone;
463 PetscInt suppSize;
464 DMPolytopeType ct;
465
466 PetscCall(DMPlexGetCellType(dm, p, &ct));
467 // Do not order tensor cells until they appear below
468 if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR || ct == DM_POLYTOPE_SEG_PRISM_TENSOR || ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) continue;
469 if (PetscBTLookupSet(bt, p)) continue;
470 PetscCall(PetscBTSet(blst, p));
471 perm[i++] = p;
472 // Check for tensor cells in the support
473 PetscCall(DMPlexGetSupport(dm, p, &supp));
474 PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
475 for (PetscInt s = 0; s < suppSize; ++s) {
476 DMPolytopeType sct;
477 PetscInt q, qq;
478
479 PetscCall(DMPlexGetCellType(dm, supp[s], &sct));
480 switch (sct) {
481 case DM_POLYTOPE_POINT_PRISM_TENSOR:
482 case DM_POLYTOPE_SEG_PRISM_TENSOR:
483 case DM_POLYTOPE_TRI_PRISM_TENSOR:
484 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
485 // If found, move up the split partner of the tensor cell, and the cell itself
486 PetscCall(DMPlexGetCone(dm, supp[s], &cone));
487 qq = supp[s];
488 q = (cone[0] == p) ? cone[1] : cone[0];
489 if (!PetscBTLookupSet(bt, q)) {
490 perm[i++] = q;
491 s = suppSize;
492 // At T-junctions, we can have an unsplit point at the other end, so also order that loop
493 {
494 const PetscInt *qsupp, *qcone;
495 PetscInt qsuppSize;
496
497 PetscCall(DMPlexGetSupport(dm, q, &qsupp));
498 PetscCall(DMPlexGetSupportSize(dm, q, &qsuppSize));
499 for (PetscInt qs = 0; qs < qsuppSize; ++qs) {
500 DMPolytopeType qsct;
501
502 PetscCall(DMPlexGetCellType(dm, qsupp[qs], &qsct));
503 switch (qsct) {
504 case DM_POLYTOPE_POINT_PRISM_TENSOR:
505 case DM_POLYTOPE_SEG_PRISM_TENSOR:
506 case DM_POLYTOPE_TRI_PRISM_TENSOR:
507 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
508 PetscCall(DMPlexGetCone(dm, qsupp[qs], &qcone));
509 if (qcone[0] == qcone[1]) {
510 if (!PetscBTLookupSet(bt, qsupp[qs])) {
511 perm[i++] = qsupp[qs];
512 qs = qsuppSize;
513 }
514 }
515 break;
516 default:
517 break;
518 }
519 }
520 }
521 }
522 if (!PetscBTLookupSet(bt, qq)) {
523 perm[i++] = qq;
524 s = suppSize;
525 }
526 break;
527 default:
528 break;
529 }
530 }
531 }
532 if (PetscDefined(USE_DEBUG)) {
533 for (PetscInt p = pStart; p < pEnd; ++p) PetscCheck(PetscBTLookup(bt, p), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " missed in permutation of [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, pStart, pEnd);
534 }
535 PetscCall(PetscBTDestroy(&bt));
536 PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart);
537 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
538 PetscCall(ISSetPermutation(permIS));
539 *permutation = permIS;
540 *blockStarts = blst;
541 PetscFunctionReturn(PETSC_SUCCESS);
542 }
543
544 // Mark the block associated with a cohesive cell p
InsertCohesiveBlock_Private(DM dm,PetscBT bt,PetscBT blst,PetscInt p,PetscInt * idx,PetscInt perm[])545 static PetscErrorCode InsertCohesiveBlock_Private(DM dm, PetscBT bt, PetscBT blst, PetscInt p, PetscInt *idx, PetscInt perm[])
546 {
547 const PetscInt *cone;
548 PetscInt cS;
549
550 PetscFunctionBegin;
551 if (PetscBTLookupSet(bt, p)) PetscFunctionReturn(PETSC_SUCCESS);
552 // Order the endcaps
553 PetscCall(DMPlexGetCone(dm, p, &cone));
554 PetscCall(DMPlexGetConeSize(dm, p, &cS));
555 if (blst) PetscCall(PetscBTSet(blst, cone[0]));
556 if (!PetscBTLookupSet(bt, cone[0])) perm[(*idx)++] = cone[0];
557 if (!PetscBTLookupSet(bt, cone[1])) perm[(*idx)++] = cone[1];
558 // Order sides
559 for (PetscInt c = 2; c < cS; ++c) PetscCall(InsertCohesiveBlock_Private(dm, bt, NULL, cone[c], idx, perm));
560 // Order cell
561 perm[(*idx)++] = p;
562 PetscFunctionReturn(PETSC_SUCCESS);
563 }
564
565 // Reorder to group split nodes
DMCreateSectionPermutation_Plex_Cohesive(DM dm,IS * permutation,PetscBT * blockStarts)566 static PetscErrorCode DMCreateSectionPermutation_Plex_Cohesive(DM dm, IS *permutation, PetscBT *blockStarts)
567 {
568 IS permIS;
569 PetscBT bt, blst;
570 PetscInt *perm;
571 PetscInt dim, pStart, pEnd, i = 0;
572
573 PetscFunctionBegin;
574 PetscCall(DMGetDimension(dm, &dim));
575 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
576 PetscCall(PetscMalloc1(pEnd - pStart, &perm));
577 PetscCall(PetscBTCreate(pEnd - pStart, &bt));
578 PetscCall(PetscBTCreate(pEnd - pStart, &blst));
579 // Add cohesive blocks
580 for (PetscInt p = pStart; p < pEnd; ++p) {
581 DMPolytopeType ct;
582
583 PetscCall(DMPlexGetCellType(dm, p, &ct));
584 switch (dim) {
585 case 2:
586 if (ct == DM_POLYTOPE_SEG_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
587 break;
588 case 3:
589 if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(InsertCohesiveBlock_Private(dm, bt, blst, p, &i, perm));
590 break;
591 default:
592 break;
593 }
594 }
595 // Add normal blocks
596 for (PetscInt p = pStart; p < pEnd; ++p) {
597 if (PetscBTLookupSet(bt, p)) continue;
598 PetscCall(PetscBTSet(blst, p));
599 perm[i++] = p;
600 }
601 if (PetscDefined(USE_DEBUG)) {
602 for (PetscInt p = pStart; p < pEnd; ++p) PetscCheck(PetscBTLookup(bt, p), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " missed in permutation of [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, pStart, pEnd);
603 }
604 PetscCall(PetscBTDestroy(&bt));
605 PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart);
606 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS));
607 PetscCall(ISSetPermutation(permIS));
608 *permutation = permIS;
609 *blockStarts = blst;
610 PetscFunctionReturn(PETSC_SUCCESS);
611 }
612
DMCreateSectionPermutation_Plex(DM dm,IS * perm,PetscBT * blockStarts)613 PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *perm, PetscBT *blockStarts)
614 {
615 DMReorderDefaultFlag reorder;
616 MatOrderingType otype;
617 PetscBool iscohesive, iscohesiveOld, isreverse;
618
619 PetscFunctionBegin;
620 PetscCall(DMReorderSectionGetDefault(dm, &reorder));
621 if (reorder != DM_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS);
622 PetscCall(DMReorderSectionGetType(dm, &otype));
623 if (!otype) PetscFunctionReturn(PETSC_SUCCESS);
624 PetscCall(PetscStrncmp(otype, "cohesive_old", 1024, &iscohesiveOld));
625 PetscCall(PetscStrncmp(otype, "cohesive", 1024, &iscohesive));
626 PetscCall(PetscStrncmp(otype, "reverse", 1024, &isreverse));
627 if (iscohesive) PetscCall(DMCreateSectionPermutation_Plex_Cohesive(dm, perm, blockStarts));
628 else if (iscohesiveOld) PetscCall(DMCreateSectionPermutation_Plex_Cohesive_Old(dm, perm, blockStarts));
629 else if (isreverse) PetscCall(DMCreateSectionPermutation_Plex_Reverse(dm, perm, blockStarts));
630 PetscFunctionReturn(PETSC_SUCCESS);
631 }
632
DMReorderSectionSetDefault_Plex(DM dm,DMReorderDefaultFlag reorder)633 PetscErrorCode DMReorderSectionSetDefault_Plex(DM dm, DMReorderDefaultFlag reorder)
634 {
635 PetscFunctionBegin;
636 dm->reorderSection = reorder;
637 PetscFunctionReturn(PETSC_SUCCESS);
638 }
639
DMReorderSectionGetDefault_Plex(DM dm,DMReorderDefaultFlag * reorder)640 PetscErrorCode DMReorderSectionGetDefault_Plex(DM dm, DMReorderDefaultFlag *reorder)
641 {
642 PetscFunctionBegin;
643 *reorder = dm->reorderSection;
644 PetscFunctionReturn(PETSC_SUCCESS);
645 }
646
DMReorderSectionSetType_Plex(DM dm,MatOrderingType reorder)647 PetscErrorCode DMReorderSectionSetType_Plex(DM dm, MatOrderingType reorder)
648 {
649 PetscFunctionBegin;
650 PetscCall(PetscFree(dm->reorderSectionType));
651 PetscCall(PetscStrallocpy(reorder, (char **)&dm->reorderSectionType));
652 PetscFunctionReturn(PETSC_SUCCESS);
653 }
654
DMReorderSectionGetType_Plex(DM dm,MatOrderingType * reorder)655 PetscErrorCode DMReorderSectionGetType_Plex(DM dm, MatOrderingType *reorder)
656 {
657 PetscFunctionBegin;
658 *reorder = dm->reorderSectionType;
659 PetscFunctionReturn(PETSC_SUCCESS);
660 }
661