1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/hashmapi.h>
3 #include <petsc/private/hashmapij.h>
4
5 const char * const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL};
6
7 /* HashIJKL */
8
9 #include <petsc/private/hashmap.h>
10
11 typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey;
12
13 #define PetscHashIJKLKeyHash(key) \
14 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \
15 PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l)))
16
17 #define PetscHashIJKLKeyEqual(k1,k2) \
18 (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0)
19
20 PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1)
21
22 static PetscSFNode _PetscInvalidSFNode = {-1, -1};
23
24 typedef struct _PetscHashIJKLRemoteKey { PetscSFNode i, j, k, l; } PetscHashIJKLRemoteKey;
25
26 #define PetscHashIJKLRemoteKeyHash(key) \
27 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index),PetscHashInt((key).j.rank + (key).j.index)), \
28 PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index),PetscHashInt((key).l.rank + (key).l.index)))
29
30 #define PetscHashIJKLRemoteKeyEqual(k1,k2) \
31 (((k1).i.rank==(k2).i.rank) ? ((k1).i.index==(k2).i.index) ? ((k1).j.rank==(k2).j.rank) ? ((k1).j.index==(k2).j.index) ? ((k1).k.rank==(k2).k.rank) ? ((k1).k.index==(k2).k.index) ? ((k1).l.rank==(k2).l.rank) ? ((k1).l.index==(k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0)
32
PETSC_HASH_MAP(HashIJKLRemote,PetscHashIJKLRemoteKey,PetscSFNode,PetscHashIJKLRemoteKeyHash,PetscHashIJKLRemoteKeyEqual,_PetscInvalidSFNode)33 PETSC_HASH_MAP(HashIJKLRemote, PetscHashIJKLRemoteKey, PetscSFNode, PetscHashIJKLRemoteKeyHash, PetscHashIJKLRemoteKeyEqual, _PetscInvalidSFNode)
34
35 static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[])
36 {
37 PetscInt i;
38
39 PetscFunctionBegin;
40 for (i = 1; i < n; ++i) {
41 PetscSFNode x = A[i];
42 PetscInt j;
43
44 for (j = i-1; j >= 0; --j) {
45 if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break;
46 A[j+1] = A[j];
47 }
48 A[j+1] = x;
49 }
50 PetscFunctionReturn(0);
51 }
52
53 /*
54 DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
55 */
DMPlexGetRawFaces_Internal(DM dm,DMPolytopeType ct,const PetscInt cone[],PetscInt * numFaces,const DMPolytopeType * faceTypes[],const PetscInt * faceSizes[],const PetscInt * faces[])56 PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
57 {
58 DMPolytopeType *typesTmp;
59 PetscInt *sizesTmp, *facesTmp;
60 PetscInt maxConeSize, maxSupportSize;
61 PetscErrorCode ierr;
62
63 PetscFunctionBegin;
64 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
65 if (cone) PetscValidIntPointer(cone, 3);
66 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
67 if (faceTypes) {ierr = DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &typesTmp);CHKERRQ(ierr);}
68 if (faceSizes) {ierr = DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &sizesTmp);CHKERRQ(ierr);}
69 if (faces) {ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);CHKERRQ(ierr);}
70 switch (ct) {
71 case DM_POLYTOPE_POINT:
72 if (numFaces) *numFaces = 0;
73 break;
74 case DM_POLYTOPE_SEGMENT:
75 if (numFaces) *numFaces = 2;
76 if (faceTypes) {
77 typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
78 *faceTypes = typesTmp;
79 }
80 if (faceSizes) {
81 sizesTmp[0] = 1; sizesTmp[1] = 1;
82 *faceSizes = sizesTmp;
83 }
84 if (faces) {
85 facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
86 *faces = facesTmp;
87 }
88 break;
89 case DM_POLYTOPE_POINT_PRISM_TENSOR:
90 if (numFaces) *numFaces = 2;
91 if (faceTypes) {
92 typesTmp[0] = DM_POLYTOPE_POINT; typesTmp[1] = DM_POLYTOPE_POINT;
93 *faceTypes = typesTmp;
94 }
95 if (faceSizes) {
96 sizesTmp[0] = 1; sizesTmp[1] = 1;
97 *faceSizes = sizesTmp;
98 }
99 if (faces) {
100 facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
101 *faces = facesTmp;
102 }
103 break;
104 case DM_POLYTOPE_TRIANGLE:
105 if (numFaces) *numFaces = 3;
106 if (faceTypes) {
107 typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT;
108 *faceTypes = typesTmp;
109 }
110 if (faceSizes) {
111 sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2;
112 *faceSizes = sizesTmp;
113 }
114 if (faces) {
115 facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
116 facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
117 facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
118 *faces = facesTmp;
119 }
120 break;
121 case DM_POLYTOPE_QUADRILATERAL:
122 /* Vertices follow right hand rule */
123 if (numFaces) *numFaces = 4;
124 if (faceTypes) {
125 typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_SEGMENT; typesTmp[3] = DM_POLYTOPE_SEGMENT;
126 *faceTypes = typesTmp;
127 }
128 if (faceSizes) {
129 sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
130 *faceSizes = sizesTmp;
131 }
132 if (faces) {
133 facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
134 facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
135 facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
136 facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
137 *faces = facesTmp;
138 }
139 break;
140 case DM_POLYTOPE_SEG_PRISM_TENSOR:
141 if (numFaces) *numFaces = 4;
142 if (faceTypes) {
143 typesTmp[0] = DM_POLYTOPE_SEGMENT; typesTmp[1] = DM_POLYTOPE_SEGMENT; typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR;
144 *faceTypes = typesTmp;
145 }
146 if (faceSizes) {
147 sizesTmp[0] = 2; sizesTmp[1] = 2; sizesTmp[2] = 2; sizesTmp[3] = 2;
148 *faceSizes = sizesTmp;
149 }
150 if (faces) {
151 facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
152 facesTmp[2] = cone[2]; facesTmp[3] = cone[3];
153 facesTmp[4] = cone[0]; facesTmp[5] = cone[2];
154 facesTmp[6] = cone[1]; facesTmp[7] = cone[3];
155 *faces = facesTmp;
156 }
157 break;
158 case DM_POLYTOPE_TETRAHEDRON:
159 /* Vertices of first face follow right hand rule and normal points away from last vertex */
160 if (numFaces) *numFaces = 4;
161 if (faceTypes) {
162 typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE; typesTmp[2] = DM_POLYTOPE_TRIANGLE; typesTmp[3] = DM_POLYTOPE_TRIANGLE;
163 *faceTypes = typesTmp;
164 }
165 if (faceSizes) {
166 sizesTmp[0] = 3; sizesTmp[1] = 3; sizesTmp[2] = 3; sizesTmp[3] = 3;
167 *faceSizes = sizesTmp;
168 }
169 if (faces) {
170 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2];
171 facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1];
172 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3];
173 facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
174 *faces = facesTmp;
175 }
176 break;
177 case DM_POLYTOPE_HEXAHEDRON:
178 /* 7--------6
179 /| /|
180 / | / |
181 4--------5 |
182 | | | |
183 | | | |
184 | 1--------2
185 | / | /
186 |/ |/
187 0--------3
188 */
189 if (numFaces) *numFaces = 6;
190 if (faceTypes) {
191 typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; typesTmp[2] = DM_POLYTOPE_QUADRILATERAL;
192 typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL; typesTmp[5] = DM_POLYTOPE_QUADRILATERAL;
193 *faceTypes = typesTmp;
194 }
195 if (faceSizes) {
196 sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
197 *faceSizes = sizesTmp;
198 }
199 if (faces) {
200 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */
201 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */
202 facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
203 facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
204 facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
205 facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
206 *faces = facesTmp;
207 }
208 break;
209 case DM_POLYTOPE_TRI_PRISM:
210 if (numFaces) *numFaces = 5;
211 if (faceTypes) {
212 typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
213 typesTmp[2] = DM_POLYTOPE_QUADRILATERAL; typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; typesTmp[4] = DM_POLYTOPE_QUADRILATERAL;
214 *faceTypes = typesTmp;
215 }
216 if (faceSizes) {
217 sizesTmp[0] = 3; sizesTmp[1] = 3;
218 sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
219 *faceSizes = sizesTmp;
220 }
221 if (faces) {
222 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; /* Bottom */
223 facesTmp[3] = cone[3]; facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; /* Top */
224 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[4]; facesTmp[9] = cone[3]; /* Back left */
225 facesTmp[10] = cone[2]; facesTmp[11] = cone[1]; facesTmp[12] = cone[5]; facesTmp[13] = cone[4]; /* Front */
226 facesTmp[14] = cone[1]; facesTmp[15] = cone[0]; facesTmp[16] = cone[3]; facesTmp[17] = cone[5]; /* Back right */
227 *faces = facesTmp;
228 }
229 break;
230 case DM_POLYTOPE_TRI_PRISM_TENSOR:
231 if (numFaces) *numFaces = 5;
232 if (faceTypes) {
233 typesTmp[0] = DM_POLYTOPE_TRIANGLE; typesTmp[1] = DM_POLYTOPE_TRIANGLE;
234 typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR;
235 *faceTypes = typesTmp;
236 }
237 if (faceSizes) {
238 sizesTmp[0] = 3; sizesTmp[1] = 3;
239 sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4;
240 *faceSizes = sizesTmp;
241 }
242 if (faces) {
243 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; /* Bottom */
244 facesTmp[3] = cone[3]; facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; /* Top */
245 facesTmp[6] = cone[0]; facesTmp[7] = cone[1]; facesTmp[8] = cone[3]; facesTmp[9] = cone[4]; /* Back left */
246 facesTmp[10] = cone[1]; facesTmp[11] = cone[2]; facesTmp[12] = cone[4]; facesTmp[13] = cone[5]; /* Back right */
247 facesTmp[14] = cone[2]; facesTmp[15] = cone[0]; facesTmp[16] = cone[5]; facesTmp[17] = cone[3]; /* Front */
248 *faces = facesTmp;
249 }
250 break;
251 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
252 /* 7--------6
253 /| /|
254 / | / |
255 4--------5 |
256 | | | |
257 | | | |
258 | 3--------2
259 | / | /
260 |/ |/
261 0--------1
262 */
263 if (numFaces) *numFaces = 6;
264 if (faceTypes) {
265 typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR;
266 typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR; typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR;
267 *faceTypes = typesTmp;
268 }
269 if (faceSizes) {
270 sizesTmp[0] = 4; sizesTmp[1] = 4; sizesTmp[2] = 4; sizesTmp[3] = 4; sizesTmp[4] = 4; sizesTmp[5] = 4;
271 *faceSizes = sizesTmp;
272 }
273 if (faces) {
274 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */
275 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */
276 facesTmp[8] = cone[0]; facesTmp[9] = cone[1]; facesTmp[10] = cone[4]; facesTmp[11] = cone[5]; /* Front */
277 facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[5]; facesTmp[15] = cone[6]; /* Right */
278 facesTmp[16] = cone[2]; facesTmp[17] = cone[3]; facesTmp[18] = cone[6]; facesTmp[19] = cone[7]; /* Back */
279 facesTmp[20] = cone[3]; facesTmp[21] = cone[0]; facesTmp[22] = cone[7]; facesTmp[23] = cone[4]; /* Left */
280 *faces = facesTmp;
281 }
282 break;
283 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]);
284 }
285 PetscFunctionReturn(0);
286 }
287
DMPlexRestoreRawFaces_Internal(DM dm,DMPolytopeType ct,const PetscInt cone[],PetscInt * numFaces,const DMPolytopeType * faceTypes[],const PetscInt * faceSizes[],const PetscInt * faces[])288 PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[])
289 {
290 PetscErrorCode ierr;
291
292 PetscFunctionBegin;
293 if (faceTypes) {ierr = DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceTypes);CHKERRQ(ierr);}
294 if (faceSizes) {ierr = DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceSizes);CHKERRQ(ierr);}
295 if (faces) {ierr = DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);CHKERRQ(ierr);}
296 PetscFunctionReturn(0);
297 }
298
299 /* This interpolates faces for cells at some stratum */
DMPlexInterpolateFaces_Internal(DM dm,PetscInt cellDepth,DM idm)300 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
301 {
302 DMLabel ctLabel;
303 PetscHashIJKL faceTable;
304 PetscInt faceTypeNum[DM_NUM_POLYTOPES];
305 PetscInt depth, d, Np, cStart, cEnd, c, fStart, fEnd;
306 PetscErrorCode ierr;
307
308 PetscFunctionBegin;
309 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
310 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr);
311 ierr = PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES);CHKERRQ(ierr);
312 ierr = DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd);CHKERRQ(ierr);
313 /* Number new faces and save face vertices in hash table */
314 ierr = DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart);CHKERRQ(ierr);
315 fEnd = fStart;
316 for (c = cStart; c < cEnd; ++c) {
317 const PetscInt *cone, *faceSizes, *faces;
318 const DMPolytopeType *faceTypes;
319 DMPolytopeType ct;
320 PetscInt numFaces, cf, foff = 0;
321
322 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr);
323 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
324 ierr = DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
325 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
326 const PetscInt faceSize = faceSizes[cf];
327 const DMPolytopeType faceType = faceTypes[cf];
328 const PetscInt *face = &faces[foff];
329 PetscHashIJKLKey key;
330 PetscHashIter iter;
331 PetscBool missing;
332
333 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
334 key.i = face[0];
335 key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
336 key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
337 key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
338 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr);
339 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
340 if (missing) {
341 ierr = PetscHashIJKLIterSet(faceTable, iter, fEnd++);CHKERRQ(ierr);
342 ++faceTypeNum[faceType];
343 }
344 }
345 ierr = DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
346 }
347 /* We need to number faces contiguously among types */
348 {
349 PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0;
350
351 for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) {if (faceTypeNum[ct]) ++numFT; faceTypeStart[ct] = 0;}
352 if (numFT > 1) {
353 ierr = PetscHashIJKLClear(faceTable);CHKERRQ(ierr);
354 faceTypeStart[0] = fStart;
355 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct-1] + faceTypeNum[ct-1];
356 for (c = cStart; c < cEnd; ++c) {
357 const PetscInt *cone, *faceSizes, *faces;
358 const DMPolytopeType *faceTypes;
359 DMPolytopeType ct;
360 PetscInt numFaces, cf, foff = 0;
361
362 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr);
363 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
364 ierr = DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
365 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
366 const PetscInt faceSize = faceSizes[cf];
367 const DMPolytopeType faceType = faceTypes[cf];
368 const PetscInt *face = &faces[foff];
369 PetscHashIJKLKey key;
370 PetscHashIter iter;
371 PetscBool missing;
372
373 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
374 key.i = face[0];
375 key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
376 key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
377 key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
378 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr);
379 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
380 if (missing) {ierr = PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++);CHKERRQ(ierr);}
381 }
382 ierr = DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
383 }
384 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) {
385 if (faceTypeStart[ct] != faceTypeStart[ct-1] + faceTypeNum[ct]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %D != %D + %D", DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct-1], faceTypeNum[ct]);
386 }
387 }
388 }
389 /* Add new points, always at the end of the numbering */
390 ierr = DMPlexGetChart(dm, NULL, &Np);CHKERRQ(ierr);
391 ierr = DMPlexSetChart(idm, 0, Np + (fEnd - fStart));CHKERRQ(ierr);
392 /* Set cone sizes */
393 /* Must create the celltype label here so that we do not automatically try to compute the types */
394 ierr = DMCreateLabel(idm, "celltype");CHKERRQ(ierr);
395 ierr = DMPlexGetCellTypeLabel(idm, &ctLabel);CHKERRQ(ierr);
396 for (d = 0; d <= depth; ++d) {
397 DMPolytopeType ct;
398 PetscInt coneSize, pStart, pEnd, p;
399
400 if (d == cellDepth) continue;
401 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr);
402 for (p = pStart; p < pEnd; ++p) {
403 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
404 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr);
405 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
406 ierr = DMPlexSetCellType(idm, p, ct);CHKERRQ(ierr);
407 }
408 }
409 for (c = cStart; c < cEnd; ++c) {
410 const PetscInt *cone, *faceSizes, *faces;
411 const DMPolytopeType *faceTypes;
412 DMPolytopeType ct;
413 PetscInt numFaces, cf, foff = 0;
414
415 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr);
416 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
417 ierr = DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
418 ierr = DMPlexSetCellType(idm, c, ct);CHKERRQ(ierr);
419 ierr = DMPlexSetConeSize(idm, c, numFaces);CHKERRQ(ierr);
420 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
421 const PetscInt faceSize = faceSizes[cf];
422 const DMPolytopeType faceType = faceTypes[cf];
423 const PetscInt *face = &faces[foff];
424 PetscHashIJKLKey key;
425 PetscHashIter iter;
426 PetscBool missing;
427 PetscInt f;
428
429 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
430 key.i = face[0];
431 key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
432 key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
433 key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
434 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr);
435 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
436 if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing face (cell %D, lf %D)", c, cf);
437 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr);
438 ierr = DMPlexSetConeSize(idm, f, faceSize);CHKERRQ(ierr);
439 ierr = DMPlexSetCellType(idm, f, faceType);CHKERRQ(ierr);
440 }
441 ierr = DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
442 }
443 ierr = DMSetUp(idm);CHKERRQ(ierr);
444 /* Initialize cones so we do not need the bash table to tell us that a cone has been set */
445 {
446 PetscSection cs;
447 PetscInt *cones, csize;
448
449 ierr = DMPlexGetConeSection(idm, &cs);CHKERRQ(ierr);
450 ierr = DMPlexGetCones(idm, &cones);CHKERRQ(ierr);
451 ierr = PetscSectionGetStorageSize(cs, &csize);CHKERRQ(ierr);
452 for (c = 0; c < csize; ++c) cones[c] = -1;
453 }
454 /* Set cones */
455 for (d = 0; d <= depth; ++d) {
456 const PetscInt *cone;
457 PetscInt pStart, pEnd, p;
458
459 if (d == cellDepth) continue;
460 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr);
461 for (p = pStart; p < pEnd; ++p) {
462 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
463 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr);
464 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr);
465 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr);
466 }
467 }
468 for (c = cStart; c < cEnd; ++c) {
469 const PetscInt *cone, *faceSizes, *faces;
470 const DMPolytopeType *faceTypes;
471 DMPolytopeType ct;
472 PetscInt numFaces, cf, foff = 0;
473
474 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr);
475 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
476 ierr = DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
477 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) {
478 DMPolytopeType faceType = faceTypes[cf];
479 const PetscInt faceSize = faceSizes[cf];
480 const PetscInt *face = &faces[foff];
481 const PetscInt *fcone;
482 PetscHashIJKLKey key;
483 PetscHashIter iter;
484 PetscBool missing;
485 PetscInt f;
486
487 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 4", faceSize);
488 key.i = face[0];
489 key.j = faceSize > 1 ? face[1] : PETSC_MAX_INT;
490 key.k = faceSize > 2 ? face[2] : PETSC_MAX_INT;
491 key.l = faceSize > 3 ? face[3] : PETSC_MAX_INT;
492 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr);
493 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
494 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr);
495 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr);
496 ierr = DMPlexGetCone(idm, f, &fcone);CHKERRQ(ierr);
497 if (fcone[0] < 0) {ierr = DMPlexSetCone(idm, f, face);CHKERRQ(ierr);}
498 /* TODO This should be unnecessary since we have autoamtic orientation */
499 {
500 /* when matching hybrid faces in 3D, only few cases are possible.
501 Face traversal however can no longer follow the usual convention, this seems a serious issue to me */
502 PetscInt tquad_map[4][4] = { {PETSC_MIN_INT, 0,PETSC_MIN_INT,PETSC_MIN_INT},
503 { -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT},
504 {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT, 1},
505 {PETSC_MIN_INT,PETSC_MIN_INT, -2,PETSC_MIN_INT} };
506 PetscInt i, i2, j;
507 const PetscInt *cone;
508 PetscInt coneSize, ornt;
509
510 /* Orient face: Do not allow reverse orientation at the first vertex */
511 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr);
512 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr);
513 if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize);
514 /* - First find the initial vertex */
515 for (i = 0; i < faceSize; ++i) if (face[0] == cone[i]) break;
516 /* If we want to compare tensor faces to regular faces, we have to flip them and take the else branch here */
517 if (faceType == DM_POLYTOPE_SEG_PRISM_TENSOR) {
518 /* find the second vertex */
519 for (i2 = 0; i2 < faceSize; ++i2) if (face[1] == cone[i2]) break;
520 switch (faceSize) {
521 case 2:
522 ornt = i ? -2 : 0;
523 break;
524 case 4:
525 ornt = tquad_map[i][i2];
526 break;
527 default: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSize, f, c);
528 }
529 } else {
530 /* Try forward comparison */
531 for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+j)%faceSize]) break;
532 if (j == faceSize) {
533 if ((faceSize == 2) && (i == 1)) ornt = -2;
534 else ornt = i;
535 } else {
536 /* Try backward comparison */
537 for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+faceSize-j)%faceSize]) break;
538 if (j == faceSize) {
539 if (i == 0) ornt = -faceSize;
540 else ornt = -i;
541 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c);
542 }
543 }
544 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr);
545 }
546 }
547 ierr = DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr);
548 }
549 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr);
550 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr);
551 ierr = DMPlexStratify(idm);CHKERRQ(ierr);
552 PetscFunctionReturn(0);
553 }
554
SortRmineRremoteByRemote_Private(PetscSF sf,PetscInt * rmine1[],PetscInt * rremote1[])555 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
556 {
557 PetscInt nleaves;
558 PetscInt nranks;
559 const PetscMPIInt *ranks=NULL;
560 const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL;
561 PetscInt n, o, r;
562 PetscErrorCode ierr;
563
564 PetscFunctionBegin;
565 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr);
566 nleaves = roffset[nranks];
567 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr);
568 for (r=0; r<nranks; r++) {
569 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
570 - to unify order with the other side */
571 o = roffset[r];
572 n = roffset[r+1] - o;
573 ierr = PetscArraycpy(&(*rmine1)[o], &rmine[o], n);CHKERRQ(ierr);
574 ierr = PetscArraycpy(&(*rremote1)[o], &rremote[o], n);CHKERRQ(ierr);
575 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr);
576 }
577 PetscFunctionReturn(0);
578 }
579
DMPlexOrientInterface_Internal(DM dm)580 PetscErrorCode DMPlexOrientInterface_Internal(DM dm)
581 {
582 /* Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
583 PetscInt masterCone[2];
584 PetscInt (*roots)[2], (*leaves)[2];
585 PetscMPIInt (*rootsRanks)[2], (*leavesRanks)[2];
586
587 PetscSF sf=NULL;
588 const PetscInt *locals=NULL;
589 const PetscSFNode *remotes=NULL;
590 PetscInt nroots, nleaves, p, c;
591 PetscInt nranks, n, o, r;
592 const PetscMPIInt *ranks=NULL;
593 const PetscInt *roffset=NULL;
594 PetscInt *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */
595 const PetscInt *cone=NULL;
596 PetscInt coneSize, ind0;
597 MPI_Comm comm;
598 PetscMPIInt rank,size;
599 PetscInt debug = 0;
600 PetscErrorCode ierr;
601
602 PetscFunctionBegin;
603 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
604 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);CHKERRQ(ierr);
605 if (nroots < 0) PetscFunctionReturn(0);
606 ierr = PetscSFSetUp(sf);CHKERRQ(ierr);
607 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr);
608 ierr = DMViewFromOptions(dm, NULL, "-before_fix_dm_view");CHKERRQ(ierr);
609 if (PetscDefined(USE_DEBUG)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);}
610 ierr = SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr);
611 ierr = PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);CHKERRQ(ierr);
612 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
613 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
614 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
615 for (p = 0; p < nroots; ++p) {
616 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
617 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
618 if (coneSize < 2) {
619 for (c = 0; c < 2; c++) {
620 roots[p][c] = -1;
621 rootsRanks[p][c] = -1;
622 }
623 continue;
624 }
625 /* Translate all points to root numbering */
626 for (c = 0; c < 2; c++) {
627 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr);
628 if (ind0 < 0) {
629 roots[p][c] = cone[c];
630 rootsRanks[p][c] = rank;
631 } else {
632 roots[p][c] = remotes[ind0].index;
633 rootsRanks[p][c] = remotes[ind0].rank;
634 }
635 }
636 }
637 for (p = 0; p < nroots; ++p) {
638 for (c = 0; c < 2; c++) {
639 leaves[p][c] = -2;
640 leavesRanks[p][c] = -2;
641 }
642 }
643 ierr = PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr);
644 ierr = PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr);
645 ierr = PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr);
646 ierr = PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr);
647 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);}
648 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Referenced roots\n");CHKERRQ(ierr);}
649 for (p = 0; p < nroots; ++p) {
650 if (leaves[p][0] < 0) continue;
651 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
652 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
653 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] %4D: cone=[%4D %4D] roots=[(%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D)]", rank, p, cone[0], cone[1], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1]);CHKERRQ(ierr);}
654 if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][1] != rootsRanks[p][1])) {
655 /* Translate these two leaves to my cone points; masterCone means desired order p's cone points */
656 for (c = 0; c < 2; c++) {
657 if (leavesRanks[p][c] == rank) {
658 /* A local leave is just taken as it is */
659 masterCone[c] = leaves[p][c];
660 continue;
661 }
662 /* Find index of rank leavesRanks[p][c] among remote ranks */
663 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */
664 ierr = PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);CHKERRQ(ierr);
665 if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): leave rank not found among remote ranks",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]);
666 if (PetscUnlikely(ranks[r] < 0 || ranks[r] >= size)) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%D c=%D commsize=%d: ranks[%D] = %d makes no sense",p,c,size,r,ranks[r]);
667 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */
668 o = roffset[r];
669 n = roffset[r+1] - o;
670 ierr = PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);CHKERRQ(ierr);
671 if (PetscUnlikely(ind0 < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): corresponding remote point not found - it seems there is missing connection in point SF!",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]);
672 /* Get the corresponding local point */
673 masterCone[c] = rmine1[o+ind0];CHKERRQ(ierr);
674 }
675 if (debug) {ierr = PetscSynchronizedPrintf(comm, " masterCone=[%4D %4D]\n", masterCone[0], masterCone[1]);CHKERRQ(ierr);}
676 /* Set the desired order of p's cone points and fix orientations accordingly */
677 /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */
678 ierr = DMPlexOrientCell(dm, p, 2, masterCone);CHKERRQ(ierr);
679 } else if (debug) {ierr = PetscSynchronizedPrintf(comm, " ==\n");CHKERRQ(ierr);}
680 }
681 if (debug) {
682 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);
683 ierr = MPI_Barrier(comm);CHKERRQ(ierr);
684 }
685 ierr = DMViewFromOptions(dm, NULL, "-after_fix_dm_view");CHKERRQ(ierr);
686 ierr = PetscFree4(roots, leaves, rootsRanks, leavesRanks);CHKERRQ(ierr);
687 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr);
688 PetscFunctionReturn(0);
689 }
690
IntArrayViewFromOptions(MPI_Comm comm,const char opt[],const char name[],const char idxname[],const char valname[],PetscInt n,const PetscInt a[])691 static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[])
692 {
693 PetscInt idx;
694 PetscMPIInt rank;
695 PetscBool flg;
696 PetscErrorCode ierr;
697
698 PetscFunctionBegin;
699 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr);
700 if (!flg) PetscFunctionReturn(0);
701 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
702 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr);
703 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);CHKERRQ(ierr);}
704 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);
705 PetscFunctionReturn(0);
706 }
707
SFNodeArrayViewFromOptions(MPI_Comm comm,const char opt[],const char name[],const char idxname[],PetscInt n,const PetscSFNode a[])708 static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[])
709 {
710 PetscInt idx;
711 PetscMPIInt rank;
712 PetscBool flg;
713 PetscErrorCode ierr;
714
715 PetscFunctionBegin;
716 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr);
717 if (!flg) PetscFunctionReturn(0);
718 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
719 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr);
720 if (idxname) {
721 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]%s %D rank %D index %D\n", rank, idxname, idx, a[idx].rank, a[idx].index);CHKERRQ(ierr);}
722 } else {
723 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);CHKERRQ(ierr);}
724 }
725 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);
726 PetscFunctionReturn(0);
727 }
728
DMPlexMapToLocalPoint(DM dm,PetscHMapIJ remotehash,PetscSFNode remotePoint,PetscInt * localPoint)729 static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint)
730 {
731 PetscSF sf;
732 const PetscInt *locals;
733 PetscMPIInt rank;
734 PetscErrorCode ierr;
735
736 PetscFunctionBegin;
737 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr);
738 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
739 ierr = PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);CHKERRQ(ierr);
740 if (remotePoint.rank == rank) {
741 *localPoint = remotePoint.index;
742 } else {
743 PetscHashIJKey key;
744 PetscInt l;
745
746 key.i = remotePoint.index;
747 key.j = remotePoint.rank;
748 ierr = PetscHMapIJGet(remotehash, key, &l);CHKERRQ(ierr);
749 if (l >= 0) {
750 *localPoint = locals[l];
751 } else PetscFunctionReturn(1);
752 }
753 PetscFunctionReturn(0);
754 }
755
DMPlexMapToGlobalPoint(DM dm,PetscInt localPoint,PetscSFNode * remotePoint)756 static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint)
757 {
758 PetscSF sf;
759 const PetscInt *locals, *rootdegree;
760 const PetscSFNode *remotes;
761 PetscInt Nl, l;
762 PetscMPIInt rank;
763 PetscErrorCode ierr;
764
765 PetscFunctionBegin;
766 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr);
767 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
768 ierr = PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);CHKERRQ(ierr);
769 if (Nl < 0) goto owned;
770 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr);
771 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr);
772 if (rootdegree[localPoint]) goto owned;
773 ierr = PetscFindInt(localPoint, Nl, locals, &l);CHKERRQ(ierr);
774 if (l < 0) PetscFunctionReturn(1);
775 *remotePoint = remotes[l];
776 PetscFunctionReturn(0);
777 owned:
778 remotePoint->rank = rank;
779 remotePoint->index = localPoint;
780 PetscFunctionReturn(0);
781 }
782
783
DMPlexPointIsShared(DM dm,PetscInt p,PetscBool * isShared)784 static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared)
785 {
786 PetscSF sf;
787 const PetscInt *locals, *rootdegree;
788 PetscInt Nl, idx;
789 PetscErrorCode ierr;
790
791 PetscFunctionBegin;
792 *isShared = PETSC_FALSE;
793 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
794 ierr = PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);CHKERRQ(ierr);
795 if (Nl < 0) PetscFunctionReturn(0);
796 ierr = PetscFindInt(p, Nl, locals, &idx);CHKERRQ(ierr);
797 if (idx >= 0) {*isShared = PETSC_TRUE; PetscFunctionReturn(0);}
798 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr);
799 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr);
800 if (rootdegree[p] > 0) *isShared = PETSC_TRUE;
801 PetscFunctionReturn(0);
802 }
803
DMPlexConeIsShared(DM dm,PetscInt p,PetscBool * isShared)804 static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared)
805 {
806 const PetscInt *cone;
807 PetscInt coneSize, c;
808 PetscBool cShared = PETSC_TRUE;
809 PetscErrorCode ierr;
810
811 PetscFunctionBegin;
812 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
813 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
814 for (c = 0; c < coneSize; ++c) {
815 PetscBool pointShared;
816
817 ierr = DMPlexPointIsShared(dm, cone[c], &pointShared);CHKERRQ(ierr);
818 cShared = (PetscBool) (cShared && pointShared);
819 }
820 *isShared = coneSize ? cShared : PETSC_FALSE;
821 PetscFunctionReturn(0);
822 }
823
DMPlexGetConeMinimum(DM dm,PetscInt p,PetscSFNode * cpmin)824 static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin)
825 {
826 const PetscInt *cone;
827 PetscInt coneSize, c;
828 PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1};
829 PetscErrorCode ierr;
830
831 PetscFunctionBegin;
832 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
833 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
834 for (c = 0; c < coneSize; ++c) {
835 PetscSFNode rcp;
836
837 ierr = DMPlexMapToGlobalPoint(dm, cone[c], &rcp);
838 if (ierr) {
839 cmin = missing;
840 } else {
841 cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin;
842 }
843 }
844 *cpmin = coneSize ? cmin : missing;
845 PetscFunctionReturn(0);
846 }
847
848 /*
849 Each shared face has an entry in the candidates array:
850 (-1, coneSize-1), {(global cone point)}
851 where the set is missing the point p which we use as the key for the face
852 */
DMPlexAddSharedFace_Private(DM dm,PetscSection candidateSection,PetscSFNode candidates[],PetscHMapIJ faceHash,PetscInt p,PetscBool debug)853 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug)
854 {
855 MPI_Comm comm;
856 const PetscInt *support;
857 PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height;
858 PetscMPIInt rank;
859 PetscErrorCode ierr;
860
861 PetscFunctionBegin;
862 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
863 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
864 ierr = DMPlexGetOverlap(dm, &overlap);CHKERRQ(ierr);
865 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
866 ierr = DMPlexGetPointHeight(dm, p, &height);CHKERRQ(ierr);
867 if (!overlap && height <= cellHeight+1) {
868 /* cells can't be shared for non-overlapping meshes */
869 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p);CHKERRQ(ierr);}
870 PetscFunctionReturn(0);
871 }
872 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr);
873 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr);
874 if (candidates) {ierr = PetscSectionGetOffset(candidateSection, p, &off);CHKERRQ(ierr);}
875 for (s = 0; s < supportSize; ++s) {
876 const PetscInt face = support[s];
877 const PetscInt *cone;
878 PetscSFNode cpmin={-1,-1}, rp={-1,-1};
879 PetscInt coneSize, c, f;
880 PetscBool isShared = PETSC_FALSE;
881 PetscHashIJKey key;
882
883 /* Only add point once */
884 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face);CHKERRQ(ierr);}
885 key.i = p;
886 key.j = face;
887 ierr = PetscHMapIJGet(faceHash, key, &f);CHKERRQ(ierr);
888 if (f >= 0) continue;
889 ierr = DMPlexConeIsShared(dm, face, &isShared);CHKERRQ(ierr);
890 ierr = DMPlexGetConeMinimum(dm, face, &cpmin);CHKERRQ(ierr);
891 ierr = DMPlexMapToGlobalPoint(dm, p, &rp);CHKERRQ(ierr);
892 if (debug) {
893 ierr = PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared);CHKERRQ(ierr);
894 ierr = PetscSynchronizedPrintf(comm, "[%d] Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index);CHKERRQ(ierr);
895 }
896 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) {
897 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr);
898 if (candidates) {
899 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank);CHKERRQ(ierr);}
900 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr);
901 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr);
902 candidates[off+idx].rank = -1;
903 candidates[off+idx++].index = coneSize-1;
904 candidates[off+idx].rank = rank;
905 candidates[off+idx++].index = face;
906 for (c = 0; c < coneSize; ++c) {
907 const PetscInt cp = cone[c];
908
909 if (cp == p) continue;
910 ierr = DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);CHKERRQ(ierr);
911 if (debug) {ierr = PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);CHKERRQ(ierr);}
912 ++idx;
913 }
914 if (debug) {ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr);}
915 } else {
916 /* Add cone size to section */
917 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face);CHKERRQ(ierr);}
918 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr);
919 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr);
920 ierr = PetscSectionAddDof(candidateSection, p, coneSize+1);CHKERRQ(ierr);
921 }
922 }
923 }
924 PetscFunctionReturn(0);
925 }
926
927 /*@
928 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation
929
930 Collective on dm
931
932 Input Parameters:
933 + dm - The interpolated DM
934 - pointSF - The initial SF without interpolated points
935
936 Output Parameter:
937 . pointSF - The SF including interpolated points
938
939 Level: developer
940
941 Note: All debugging for this process can be turned on with the options: -dm_interp_pre_view -petscsf_interp_pre_view -petscsection_interp_candidate_view -petscsection_interp_candidate_remote_view -petscsection_interp_claim_view -petscsf_interp_pre_view -dmplex_interp_debug
942
943 .seealso: DMPlexInterpolate(), DMPlexUninterpolate()
944 @*/
DMPlexInterpolatePointSF(DM dm,PetscSF pointSF)945 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF)
946 {
947 MPI_Comm comm;
948 PetscHMapIJ remoteHash;
949 PetscHMapI claimshash;
950 PetscSection candidateSection, candidateRemoteSection, claimSection;
951 PetscSFNode *candidates, *candidatesRemote, *claims;
952 const PetscInt *localPoints, *rootdegree;
953 const PetscSFNode *remotePoints;
954 PetscInt ov, Nr, r, Nl, l;
955 PetscInt candidatesSize, candidatesRemoteSize, claimsSize;
956 PetscBool flg, debug = PETSC_FALSE;
957 PetscMPIInt rank;
958 PetscErrorCode ierr;
959
960 PetscFunctionBegin;
961 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
962 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 3);
963 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr);
964 if (!flg) PetscFunctionReturn(0);
965 /* Set initial SF so that lower level queries work */
966 ierr = DMSetPointSF(dm, pointSF);CHKERRQ(ierr);
967 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
968 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
969 ierr = DMPlexGetOverlap(dm, &ov);CHKERRQ(ierr);
970 if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet");
971 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr);
972 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr);
973 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr);
974 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr);
975 /* Step 0: Precalculations */
976 ierr = PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);CHKERRQ(ierr);
977 if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set");
978 ierr = PetscHMapIJCreate(&remoteHash);CHKERRQ(ierr);
979 for (l = 0; l < Nl; ++l) {
980 PetscHashIJKey key;
981 key.i = remotePoints[l].index;
982 key.j = remotePoints[l].rank;
983 ierr = PetscHMapIJSet(remoteHash, key, l);CHKERRQ(ierr);
984 }
985 /* Compute root degree to identify shared points */
986 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr);
987 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr);
988 ierr = IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);CHKERRQ(ierr);
989 /*
990 1) Loop over each leaf point $p$ at depth $d$ in the SF
991 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which
992 \begin{itemize}
993 \item all cone points of $f$ are shared
994 \item $p$ is the cone point with smallest canonical number
995 \end{itemize}
996 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$
997 \item At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared \label{alg:rootStep} and choose the root face
998 \item Send the root face from the root back to all leaf process
999 \item Leaf processes add the shared face to the SF
1000 */
1001 /* Step 1: Construct section+SFNode array
1002 The section has entries for all shared faces for which we have a leaf point in the cone
1003 The array holds candidate shared faces, each face is refered to by the leaf point */
1004 ierr = PetscSectionCreate(comm, &candidateSection);CHKERRQ(ierr);
1005 ierr = PetscSectionSetChart(candidateSection, 0, Nr);CHKERRQ(ierr);
1006 {
1007 PetscHMapIJ faceHash;
1008
1009 ierr = PetscHMapIJCreate(&faceHash);CHKERRQ(ierr);
1010 for (l = 0; l < Nl; ++l) {
1011 const PetscInt p = localPoints[l];
1012
1013 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p);CHKERRQ(ierr);}
1014 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);CHKERRQ(ierr);
1015 }
1016 ierr = PetscHMapIJClear(faceHash);CHKERRQ(ierr);
1017 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr);
1018 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr);
1019 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr);
1020 for (l = 0; l < Nl; ++l) {
1021 const PetscInt p = localPoints[l];
1022
1023 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p);CHKERRQ(ierr);}
1024 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);CHKERRQ(ierr);
1025 }
1026 ierr = PetscHMapIJDestroy(&faceHash);CHKERRQ(ierr);
1027 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);}
1028 }
1029 ierr = PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");CHKERRQ(ierr);
1030 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr);
1031 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr);
1032 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */
1033 /* Note that this section is indexed by offsets into leaves, not by point number */
1034 {
1035 PetscSF sfMulti, sfInverse, sfCandidates;
1036 PetscInt *remoteOffsets;
1037
1038 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr);
1039 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr);
1040 ierr = PetscSectionCreate(comm, &candidateRemoteSection);CHKERRQ(ierr);
1041 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);CHKERRQ(ierr);
1042 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);CHKERRQ(ierr);
1043 ierr = PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);CHKERRQ(ierr);
1044 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr);
1045 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr);
1046 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr);
1047 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr);
1048 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr);
1049 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
1050
1051 ierr = PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");CHKERRQ(ierr);
1052 ierr = PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr);
1053 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr);
1054 }
1055 /* Step 3: At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared and choose the root face */
1056 {
1057 PetscHashIJKLRemote faceTable;
1058 PetscInt idx, idx2;
1059
1060 ierr = PetscHashIJKLRemoteCreate(&faceTable);CHKERRQ(ierr);
1061 /* There is a section point for every leaf attached to a given root point */
1062 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) {
1063 PetscInt deg;
1064
1065 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) {
1066 PetscInt offset, dof, d;
1067
1068 ierr = PetscSectionGetDof(candidateRemoteSection, idx, &dof);CHKERRQ(ierr);
1069 ierr = PetscSectionGetOffset(candidateRemoteSection, idx, &offset);CHKERRQ(ierr);
1070 /* dof may include many faces from the remote process */
1071 for (d = 0; d < dof; ++d) {
1072 const PetscInt hidx = offset+d;
1073 const PetscInt Np = candidatesRemote[hidx].index+1;
1074 const PetscSFNode rface = candidatesRemote[hidx+1];
1075 const PetscSFNode *fcone = &candidatesRemote[hidx+2];
1076 PetscSFNode fcp0;
1077 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1078 const PetscInt *join = NULL;
1079 PetscHashIJKLRemoteKey key;
1080 PetscHashIter iter;
1081 PetscBool missing;
1082 PetscInt points[1024], p, joinSize;
1083
1084 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking face (%D, %D) at (%D, %D, %D) with cone size %D\n", rank, rface.rank, rface.index, r, idx, d, Np);CHKERRQ(ierr);}
1085 if (Np > 4) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%D, %D) at (%D, %D, %D) with %D cone points", rface.rank, rface.index, r, idx, d, Np);
1086 fcp0.rank = rank;
1087 fcp0.index = r;
1088 d += Np;
1089 /* Put remote face in hash table */
1090 key.i = fcp0;
1091 key.j = fcone[0];
1092 key.k = Np > 2 ? fcone[1] : pmax;
1093 key.l = Np > 3 ? fcone[2] : pmax;
1094 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr);
1095 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
1096 if (missing) {
1097 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);}
1098 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr);
1099 } else {
1100 PetscSFNode oface;
1101
1102 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr);
1103 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) {
1104 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);}
1105 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr);
1106 }
1107 }
1108 /* Check for local face */
1109 points[0] = r;
1110 for (p = 1; p < Np; ++p) {
1111 ierr = DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]);
1112 if (ierr) break; /* We got a point not in our overlap */
1113 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p]);CHKERRQ(ierr);}
1114 }
1115 if (ierr) continue;
1116 ierr = DMPlexGetJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr);
1117 if (joinSize == 1) {
1118 PetscSFNode lface;
1119 PetscSFNode oface;
1120
1121 /* Always replace with local face */
1122 lface.rank = rank;
1123 lface.index = join[0];
1124 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr);
1125 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing (%D, %D) with local face (%D, %D)\n", rank, oface.index, oface.rank, lface.index, lface.rank);CHKERRQ(ierr);}
1126 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, lface);CHKERRQ(ierr);
1127 }
1128 ierr = DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr);
1129 }
1130 }
1131 /* Put back faces for this root */
1132 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) {
1133 PetscInt offset, dof, d;
1134
1135 ierr = PetscSectionGetDof(candidateRemoteSection, idx2, &dof);CHKERRQ(ierr);
1136 ierr = PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);CHKERRQ(ierr);
1137 /* dof may include many faces from the remote process */
1138 for (d = 0; d < dof; ++d) {
1139 const PetscInt hidx = offset+d;
1140 const PetscInt Np = candidatesRemote[hidx].index+1;
1141 const PetscSFNode *fcone = &candidatesRemote[hidx+2];
1142 PetscSFNode fcp0;
1143 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT};
1144 PetscHashIJKLRemoteKey key;
1145 PetscHashIter iter;
1146 PetscBool missing;
1147
1148 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx);CHKERRQ(ierr);}
1149 if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np);
1150 fcp0.rank = rank;
1151 fcp0.index = r;
1152 d += Np;
1153 /* Find remote face in hash table */
1154 key.i = fcp0;
1155 key.j = fcone[0];
1156 key.k = Np > 2 ? fcone[1] : pmax;
1157 key.l = Np > 3 ? fcone[2] : pmax;
1158 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr);
1159 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] key (%D, %D) (%D, %D) (%D, %D) (%D, %D)\n", rank, key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index);CHKERRQ(ierr);}
1160 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr);
1161 if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2);
1162 else {ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);CHKERRQ(ierr);}
1163 }
1164 }
1165 }
1166 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);}
1167 ierr = PetscHashIJKLRemoteDestroy(&faceTable);CHKERRQ(ierr);
1168 }
1169 /* Step 4: Push back owned faces */
1170 {
1171 PetscSF sfMulti, sfClaims, sfPointNew;
1172 PetscSFNode *remotePointsNew;
1173 PetscInt *remoteOffsets, *localPointsNew;
1174 PetscInt pStart, pEnd, r, NlNew, p;
1175
1176 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */
1177 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr);
1178 ierr = PetscSectionCreate(comm, &claimSection);CHKERRQ(ierr);
1179 ierr = PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);CHKERRQ(ierr);
1180 ierr = PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr);
1181 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr);
1182 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr);
1183 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1;
1184 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr);
1185 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr);
1186 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr);
1187 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
1188 ierr = PetscObjectSetName((PetscObject) claimSection, "Claim Section");CHKERRQ(ierr);
1189 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr);
1190 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr);
1191 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */
1192 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */
1193 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr);
1194 for (r = 0; r < Nr; ++r) {
1195 PetscInt dof, off, d;
1196
1197 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r);CHKERRQ(ierr);}
1198 ierr = PetscSectionGetDof(candidateSection, r, &dof);CHKERRQ(ierr);
1199 ierr = PetscSectionGetOffset(candidateSection, r, &off);CHKERRQ(ierr);
1200 for (d = 0; d < dof;) {
1201 if (claims[off+d].rank >= 0) {
1202 const PetscInt faceInd = off+d;
1203 const PetscInt Np = candidates[off+d].index;
1204 const PetscInt *join = NULL;
1205 PetscInt joinSize, points[1024], c;
1206
1207 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);CHKERRQ(ierr);}
1208 points[0] = r;
1209 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);}
1210 for (c = 0, d += 2; c < Np; ++c, ++d) {
1211 ierr = DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);CHKERRQ(ierr);
1212 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);}
1213 }
1214 ierr = DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr);
1215 if (joinSize == 1) {
1216 if (claims[faceInd].rank == rank) {
1217 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0]);CHKERRQ(ierr);}
1218 } else {
1219 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);}
1220 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr);
1221 }
1222 } else {
1223 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank);CHKERRQ(ierr);}
1224 }
1225 ierr = DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr);
1226 } else {
1227 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r);CHKERRQ(ierr);}
1228 d += claims[off+d].index+1;
1229 }
1230 }
1231 }
1232 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);}
1233 /* Step 6) Create new pointSF from hashed claims */
1234 ierr = PetscHMapIGetSize(claimshash, &NlNew);CHKERRQ(ierr);
1235 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
1236 ierr = PetscMalloc1(Nl + NlNew, &localPointsNew);CHKERRQ(ierr);
1237 ierr = PetscMalloc1(Nl + NlNew, &remotePointsNew);CHKERRQ(ierr);
1238 for (l = 0; l < Nl; ++l) {
1239 localPointsNew[l] = localPoints[l];
1240 remotePointsNew[l].index = remotePoints[l].index;
1241 remotePointsNew[l].rank = remotePoints[l].rank;
1242 }
1243 p = Nl;
1244 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr);
1245 /* We sort new points, and assume they are numbered after all existing points */
1246 ierr = PetscSortInt(NlNew, &localPointsNew[Nl]);CHKERRQ(ierr);
1247 for (p = Nl; p < Nl + NlNew; ++p) {
1248 PetscInt off;
1249 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &off);CHKERRQ(ierr);
1250 if (claims[off].rank < 0 || claims[off].index < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %D, (%D, %D)", localPointsNew[p], claims[off].rank, claims[off].index);
1251 remotePointsNew[p] = claims[off];
1252 }
1253 ierr = PetscSFCreate(comm, &sfPointNew);CHKERRQ(ierr);
1254 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr);
1255 ierr = PetscSFSetUp(sfPointNew);CHKERRQ(ierr);
1256 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr);
1257 ierr = PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");CHKERRQ(ierr);
1258 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr);
1259 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr);
1260 }
1261 ierr = PetscHMapIJDestroy(&remoteHash);CHKERRQ(ierr);
1262 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr);
1263 ierr = PetscSectionDestroy(&candidateRemoteSection);CHKERRQ(ierr);
1264 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr);
1265 ierr = PetscFree(candidates);CHKERRQ(ierr);
1266 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr);
1267 ierr = PetscFree(claims);CHKERRQ(ierr);
1268 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr);
1269 PetscFunctionReturn(0);
1270 }
1271
1272 /*@
1273 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
1274
1275 Collective on dm
1276
1277 Input Parameters:
1278 + dm - The DMPlex object with only cells and vertices
1279 - dmInt - The interpolated DM
1280
1281 Output Parameter:
1282 . dmInt - The complete DMPlex object
1283
1284 Level: intermediate
1285
1286 Notes:
1287 It does not copy over the coordinates.
1288
1289 Developer Notes:
1290 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL.
1291
1292 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1293 @*/
DMPlexInterpolate(DM dm,DM * dmInt)1294 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
1295 {
1296 DMPlexInterpolatedFlag interpolated;
1297 DM idm, odm = dm;
1298 PetscSF sfPoint;
1299 PetscInt depth, dim, d;
1300 const char *name;
1301 PetscBool flg=PETSC_TRUE;
1302 PetscErrorCode ierr;
1303
1304 PetscFunctionBegin;
1305 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1306 PetscValidPointer(dmInt, 2);
1307 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr);
1308 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
1309 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1310 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr);
1311 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1312 if (interpolated == DMPLEX_INTERPOLATED_FULL) {
1313 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
1314 idm = dm;
1315 } else {
1316 for (d = 1; d < dim; ++d) {
1317 /* Create interpolated mesh */
1318 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr);
1319 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr);
1320 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr);
1321 if (depth > 0) {
1322 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr);
1323 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr);
1324 {
1325 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */
1326 PetscInt nroots;
1327 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
1328 if (nroots >= 0) {ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr);}
1329 }
1330 }
1331 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);}
1332 odm = idm;
1333 }
1334 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr);
1335 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr);
1336 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
1337 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr);
1338 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr);
1339 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);}
1340 }
1341 {
1342 PetscBool isper;
1343 const PetscReal *maxCell, *L;
1344 const DMBoundaryType *bd;
1345
1346 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr);
1347 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr);
1348 }
1349 /* This function makes the mesh fully interpolated on all ranks */
1350 {
1351 DM_Plex *plex = (DM_Plex *) idm->data;
1352 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL;
1353 }
1354 *dmInt = idm;
1355 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr);
1356 PetscFunctionReturn(0);
1357 }
1358
1359 /*@
1360 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
1361
1362 Collective on dmA
1363
1364 Input Parameter:
1365 . dmA - The DMPlex object with initial coordinates
1366
1367 Output Parameter:
1368 . dmB - The DMPlex object with copied coordinates
1369
1370 Level: intermediate
1371
1372 Note: This is typically used when adding pieces other than vertices to a mesh
1373
1374 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
1375 @*/
DMPlexCopyCoordinates(DM dmA,DM dmB)1376 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
1377 {
1378 Vec coordinatesA, coordinatesB;
1379 VecType vtype;
1380 PetscSection coordSectionA, coordSectionB;
1381 PetscScalar *coordsA, *coordsB;
1382 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
1383 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim;
1384 PetscBool lc = PETSC_FALSE;
1385 PetscErrorCode ierr;
1386
1387 PetscFunctionBegin;
1388 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1);
1389 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2);
1390 if (dmA == dmB) PetscFunctionReturn(0);
1391 ierr = DMGetCoordinateDim(dmA, &cdim);CHKERRQ(ierr);
1392 ierr = DMSetCoordinateDim(dmB, cdim);CHKERRQ(ierr);
1393 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr);
1394 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr);
1395 if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB);
1396 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr);
1397 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr);
1398 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr);
1399 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr);
1400 if (coordSectionA == coordSectionB) PetscFunctionReturn(0);
1401 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr);
1402 if (!Nf) PetscFunctionReturn(0);
1403 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf);
1404 if (!coordSectionB) {
1405 PetscInt dim;
1406
1407 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr);
1408 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr);
1409 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr);
1410 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr);
1411 }
1412 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr);
1413 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr);
1414 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr);
1415 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr);
1416 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */
1417 if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D in the second DM", cEndA-cStartA, cEndB-cStartB);
1418 cS = cS - cStartA + cStartB;
1419 cE = vEndB;
1420 lc = PETSC_TRUE;
1421 } else {
1422 cS = vStartB;
1423 cE = vEndB;
1424 }
1425 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr);
1426 for (v = vStartB; v < vEndB; ++v) {
1427 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr);
1428 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr);
1429 }
1430 if (lc) { /* localized coordinates */
1431 PetscInt c;
1432
1433 for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1434 PetscInt dof;
1435
1436 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr);
1437 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr);
1438 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr);
1439 }
1440 }
1441 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr);
1442 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr);
1443 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr);
1444 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr);
1445 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr);
1446 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr);
1447 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr);
1448 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr);
1449 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr);
1450 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr);
1451 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr);
1452 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr);
1453 for (v = 0; v < vEndB-vStartB; ++v) {
1454 PetscInt offA, offB;
1455
1456 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr);
1457 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr);
1458 for (d = 0; d < spaceDim; ++d) {
1459 coordsB[offB+d] = coordsA[offA+d];
1460 }
1461 }
1462 if (lc) { /* localized coordinates */
1463 PetscInt c;
1464
1465 for (c = cS-cStartB; c < cEndB-cStartB; c++) {
1466 PetscInt dof, offA, offB;
1467
1468 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr);
1469 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr);
1470 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr);
1471 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr);
1472 }
1473 }
1474 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr);
1475 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr);
1476 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr);
1477 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr);
1478 PetscFunctionReturn(0);
1479 }
1480
1481 /*@
1482 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
1483
1484 Collective on dm
1485
1486 Input Parameter:
1487 . dm - The complete DMPlex object
1488
1489 Output Parameter:
1490 . dmUnint - The DMPlex object with only cells and vertices
1491
1492 Level: intermediate
1493
1494 Notes:
1495 It does not copy over the coordinates.
1496
1497 Developer Notes:
1498 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.
1499
1500 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates()
1501 @*/
DMPlexUninterpolate(DM dm,DM * dmUnint)1502 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
1503 {
1504 DMPlexInterpolatedFlag interpolated;
1505 DM udm;
1506 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
1507 PetscErrorCode ierr;
1508
1509 PetscFunctionBegin;
1510 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1511 PetscValidPointer(dmUnint, 2);
1512 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1513 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr);
1514 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
1515 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) {
1516 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */
1517 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
1518 *dmUnint = dm;
1519 PetscFunctionReturn(0);
1520 }
1521 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
1522 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1523 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr);
1524 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr);
1525 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr);
1526 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr);
1527 for (c = cStart; c < cEnd; ++c) {
1528 PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1529
1530 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1531 for (cl = 0; cl < closureSize*2; cl += 2) {
1532 const PetscInt p = closure[cl];
1533
1534 if ((p >= vStart) && (p < vEnd)) ++coneSize;
1535 }
1536 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1537 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr);
1538 maxConeSize = PetscMax(maxConeSize, coneSize);
1539 }
1540 ierr = DMSetUp(udm);CHKERRQ(ierr);
1541 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr);
1542 for (c = cStart; c < cEnd; ++c) {
1543 PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
1544
1545 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1546 for (cl = 0; cl < closureSize*2; cl += 2) {
1547 const PetscInt p = closure[cl];
1548
1549 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
1550 }
1551 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1552 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr);
1553 }
1554 ierr = PetscFree(cone);CHKERRQ(ierr);
1555 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr);
1556 ierr = DMPlexStratify(udm);CHKERRQ(ierr);
1557 /* Reduce SF */
1558 {
1559 PetscSF sfPoint, sfPointUn;
1560 const PetscSFNode *remotePoints;
1561 const PetscInt *localPoints;
1562 PetscSFNode *remotePointsUn;
1563 PetscInt *localPointsUn;
1564 PetscInt vEnd, numRoots, numLeaves, l;
1565 PetscInt numLeavesUn = 0, n = 0;
1566 PetscErrorCode ierr;
1567
1568 /* Get original SF information */
1569 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr);
1570 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr);
1571 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr);
1572 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr);
1573 /* Allocate space for cells and vertices */
1574 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++;
1575 /* Fill in leaves */
1576 if (vEnd >= 0) {
1577 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr);
1578 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr);
1579 for (l = 0; l < numLeaves; l++) {
1580 if (localPoints[l] < vEnd) {
1581 localPointsUn[n] = localPoints[l];
1582 remotePointsUn[n].rank = remotePoints[l].rank;
1583 remotePointsUn[n].index = remotePoints[l].index;
1584 ++n;
1585 }
1586 }
1587 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn);
1588 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr);
1589 }
1590 }
1591 {
1592 PetscBool isper;
1593 const PetscReal *maxCell, *L;
1594 const DMBoundaryType *bd;
1595
1596 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr);
1597 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr);
1598 }
1599 /* This function makes the mesh fully uninterpolated on all ranks */
1600 {
1601 DM_Plex *plex = (DM_Plex *) udm->data;
1602 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE;
1603 }
1604 *dmUnint = udm;
1605 PetscFunctionReturn(0);
1606 }
1607
DMPlexIsInterpolated_Internal(DM dm,DMPlexInterpolatedFlag * interpolated)1608 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated)
1609 {
1610 PetscInt coneSize, depth, dim, h, p, pStart, pEnd;
1611 MPI_Comm comm;
1612 PetscErrorCode ierr;
1613
1614 PetscFunctionBegin;
1615 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
1616 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
1617 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1618
1619 if (depth == dim) {
1620 *interpolated = DMPLEX_INTERPOLATED_FULL;
1621 if (!dim) goto finish;
1622
1623 /* Check points at height = dim are vertices (have no cones) */
1624 ierr = DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);CHKERRQ(ierr);
1625 for (p=pStart; p<pEnd; p++) {
1626 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
1627 if (coneSize) {
1628 *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1629 goto finish;
1630 }
1631 }
1632
1633 /* Check points at height < dim have cones */
1634 for (h=0; h<dim; h++) {
1635 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr);
1636 for (p=pStart; p<pEnd; p++) {
1637 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
1638 if (!coneSize) {
1639 *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1640 goto finish;
1641 }
1642 }
1643 }
1644 } else if (depth == 1) {
1645 *interpolated = DMPLEX_INTERPOLATED_NONE;
1646 } else {
1647 *interpolated = DMPLEX_INTERPOLATED_PARTIAL;
1648 }
1649 finish:
1650 PetscFunctionReturn(0);
1651 }
1652
1653 /*@
1654 DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated.
1655
1656 Not Collective
1657
1658 Input Parameter:
1659 . dm - The DM object
1660
1661 Output Parameter:
1662 . interpolated - Flag whether the DM is interpolated
1663
1664 Level: intermediate
1665
1666 Notes:
1667 Unlike DMPlexIsInterpolatedCollective(), this is NOT collective
1668 so the results can be different on different ranks in special cases.
1669 However, DMPlexInterpolate() guarantees the result is the same on all.
1670
1671 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED.
1672
1673 Developer Notes:
1674 Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID.
1675
1676 If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called.
1677 It checks the actual topology and sets plex->interpolated on each rank separately to one of
1678 DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL.
1679
1680 If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated.
1681
1682 DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL,
1683 and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE.
1684
1685 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective()
1686 @*/
DMPlexIsInterpolated(DM dm,DMPlexInterpolatedFlag * interpolated)1687 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated)
1688 {
1689 DM_Plex *plex = (DM_Plex *) dm->data;
1690 PetscErrorCode ierr;
1691
1692 PetscFunctionBegin;
1693 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1694 PetscValidPointer(interpolated,2);
1695 if (plex->interpolated < 0) {
1696 ierr = DMPlexIsInterpolated_Internal(dm, &plex->interpolated);CHKERRQ(ierr);
1697 } else if (PetscDefined (USE_DEBUG)) {
1698 DMPlexInterpolatedFlag flg;
1699
1700 ierr = DMPlexIsInterpolated_Internal(dm, &flg);CHKERRQ(ierr);
1701 if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]);
1702 }
1703 *interpolated = plex->interpolated;
1704 PetscFunctionReturn(0);
1705 }
1706
1707 /*@
1708 DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner).
1709
1710 Collective
1711
1712 Input Parameter:
1713 . dm - The DM object
1714
1715 Output Parameter:
1716 . interpolated - Flag whether the DM is interpolated
1717
1718 Level: intermediate
1719
1720 Notes:
1721 Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks.
1722
1723 This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks.
1724
1725 Developer Notes:
1726 Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID.
1727
1728 If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated.
1729 MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned;
1730 if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED,
1731 otherwise sets plex->interpolatedCollective = plex->interpolated.
1732
1733 If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective.
1734
1735 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated()
1736 @*/
DMPlexIsInterpolatedCollective(DM dm,DMPlexInterpolatedFlag * interpolated)1737 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated)
1738 {
1739 DM_Plex *plex = (DM_Plex *) dm->data;
1740 PetscBool debug=PETSC_FALSE;
1741 PetscErrorCode ierr;
1742
1743 PetscFunctionBegin;
1744 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1745 PetscValidPointer(interpolated,2);
1746 ierr = PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);CHKERRQ(ierr);
1747 if (plex->interpolatedCollective < 0) {
1748 DMPlexInterpolatedFlag min, max;
1749 MPI_Comm comm;
1750
1751 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
1752 ierr = DMPlexIsInterpolated(dm, &plex->interpolatedCollective);CHKERRQ(ierr);
1753 ierr = MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);CHKERRQ(ierr);
1754 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRQ(ierr);
1755 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED;
1756 if (debug) {
1757 PetscMPIInt rank;
1758
1759 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1760 ierr = PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);CHKERRQ(ierr);
1761 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr);
1762 }
1763 }
1764 *interpolated = plex->interpolatedCollective;
1765 PetscFunctionReturn(0);
1766 }
1767