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