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)); 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) 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 else 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 PetscCall(PetscSFDestroy(&sfPointNew)); 1265 PetscCall(PetscHMapIDestroy(&claimshash)); 1266 } 1267 PetscCall(PetscHMapIJDestroy(&remoteHash)); 1268 PetscCall(PetscSectionDestroy(&candidateSection)); 1269 PetscCall(PetscSectionDestroy(&candidateRemoteSection)); 1270 PetscCall(PetscSectionDestroy(&claimSection)); 1271 PetscCall(PetscFree(candidates)); 1272 PetscCall(PetscFree(candidatesRemote)); 1273 PetscCall(PetscFree(claims)); 1274 PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0)); 1275 PetscFunctionReturn(0); 1276 } 1277 1278 /*@ 1279 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1280 1281 Collective on dm 1282 1283 Input Parameters: 1284 + dm - The DMPlex object with only cells and vertices 1285 - dmInt - The interpolated DM 1286 1287 Output Parameter: 1288 . dmInt - The complete DMPlex object 1289 1290 Level: intermediate 1291 1292 Notes: 1293 Labels and coordinates are copied. 1294 1295 Developer Notes: 1296 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL. 1297 1298 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1299 @*/ 1300 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1301 { 1302 DMPlexInterpolatedFlag interpolated; 1303 DM idm, odm = dm; 1304 PetscSF sfPoint; 1305 PetscInt depth, dim, d; 1306 const char *name; 1307 PetscBool flg=PETSC_TRUE; 1308 1309 PetscFunctionBegin; 1310 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1311 PetscValidPointer(dmInt, 2); 1312 PetscCall(PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0)); 1313 PetscCall(DMPlexGetDepth(dm, &depth)); 1314 PetscCall(DMGetDimension(dm, &dim)); 1315 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1316 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1317 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1318 PetscCall(PetscObjectReference((PetscObject) dm)); 1319 idm = dm; 1320 } else { 1321 for (d = 1; d < dim; ++d) { 1322 /* Create interpolated mesh */ 1323 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1324 PetscCall(DMSetType(idm, DMPLEX)); 1325 PetscCall(DMSetDimension(idm, dim)); 1326 if (depth > 0) { 1327 PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm)); 1328 PetscCall(DMGetPointSF(odm, &sfPoint)); 1329 { 1330 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1331 PetscInt nroots; 1332 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1333 if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint)); 1334 } 1335 } 1336 if (odm != dm) PetscCall(DMDestroy(&odm)); 1337 odm = idm; 1338 } 1339 PetscCall(PetscObjectGetName((PetscObject) dm, &name)); 1340 PetscCall(PetscObjectSetName((PetscObject) idm, name)); 1341 PetscCall(DMPlexCopyCoordinates(dm, idm)); 1342 PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1343 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL)); 1344 if (flg) PetscCall(DMPlexOrientInterface_Internal(idm)); 1345 } 1346 /* This function makes the mesh fully interpolated on all ranks */ 1347 { 1348 DM_Plex *plex = (DM_Plex *) idm->data; 1349 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1350 } 1351 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, idm)); 1352 *dmInt = idm; 1353 PetscCall(PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0)); 1354 PetscFunctionReturn(0); 1355 } 1356 1357 /*@ 1358 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1359 1360 Collective on dmA 1361 1362 Input Parameter: 1363 . dmA - The DMPlex object with initial coordinates 1364 1365 Output Parameter: 1366 . dmB - The DMPlex object with copied coordinates 1367 1368 Level: intermediate 1369 1370 Note: This is typically used when adding pieces other than vertices to a mesh 1371 1372 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1373 @*/ 1374 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1375 { 1376 Vec coordinatesA, coordinatesB; 1377 VecType vtype; 1378 PetscSection coordSectionA, coordSectionB; 1379 PetscScalar *coordsA, *coordsB; 1380 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1381 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1382 PetscBool lc = PETSC_FALSE; 1383 1384 PetscFunctionBegin; 1385 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1386 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1387 if (dmA == dmB) PetscFunctionReturn(0); 1388 PetscCall(DMGetCoordinateDim(dmA, &cdim)); 1389 PetscCall(DMSetCoordinateDim(dmB, cdim)); 1390 PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA)); 1391 PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB)); 1392 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); 1393 /* Copy over discretization if it exists */ 1394 { 1395 DM cdmA, cdmB; 1396 PetscDS dsA, dsB; 1397 PetscObject objA, objB; 1398 PetscClassId idA, idB; 1399 const PetscScalar *constants; 1400 PetscInt cdim, Nc; 1401 1402 PetscCall(DMGetCoordinateDM(dmA, &cdmA)); 1403 PetscCall(DMGetCoordinateDM(dmB, &cdmB)); 1404 PetscCall(DMGetField(cdmA, 0, NULL, &objA)); 1405 PetscCall(DMGetField(cdmB, 0, NULL, &objB)); 1406 PetscCall(PetscObjectGetClassId(objA, &idA)); 1407 PetscCall(PetscObjectGetClassId(objB, &idB)); 1408 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1409 PetscCall(DMSetField(cdmB, 0, NULL, objA)); 1410 PetscCall(DMCreateDS(cdmB)); 1411 PetscCall(DMGetDS(cdmA, &dsA)); 1412 PetscCall(DMGetDS(cdmB, &dsB)); 1413 PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim)); 1414 PetscCall(PetscDSSetCoordinateDimension(dsB, cdim)); 1415 PetscCall(PetscDSGetConstants(dsA, &Nc, &constants)); 1416 PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants)); 1417 } 1418 } 1419 PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA)); 1420 PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB)); 1421 PetscCall(DMGetCoordinateSection(dmA, &coordSectionA)); 1422 PetscCall(DMGetCoordinateSection(dmB, &coordSectionB)); 1423 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1424 PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf)); 1425 if (!Nf) PetscFunctionReturn(0); 1426 PetscCheck(Nf <= 1,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf); 1427 if (!coordSectionB) { 1428 PetscInt dim; 1429 1430 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB)); 1431 PetscCall(DMGetCoordinateDim(dmA, &dim)); 1432 PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB)); 1433 PetscCall(PetscObjectDereference((PetscObject) coordSectionB)); 1434 } 1435 PetscCall(PetscSectionSetNumFields(coordSectionB, 1)); 1436 PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim)); 1437 PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim)); 1438 PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE)); 1439 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1440 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); 1441 cS = cS - cStartA + cStartB; 1442 cE = vEndB; 1443 lc = PETSC_TRUE; 1444 } else { 1445 cS = vStartB; 1446 cE = vEndB; 1447 } 1448 PetscCall(PetscSectionSetChart(coordSectionB, cS, cE)); 1449 for (v = vStartB; v < vEndB; ++v) { 1450 PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim)); 1451 PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim)); 1452 } 1453 if (lc) { /* localized coordinates */ 1454 PetscInt c; 1455 1456 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1457 PetscInt dof; 1458 1459 PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1460 PetscCall(PetscSectionSetDof(coordSectionB, c + cStartB, dof)); 1461 PetscCall(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof)); 1462 } 1463 } 1464 PetscCall(PetscSectionSetUp(coordSectionB)); 1465 PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB)); 1466 PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA)); 1467 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB)); 1468 PetscCall(PetscObjectSetName((PetscObject) coordinatesB, "coordinates")); 1469 PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE)); 1470 PetscCall(VecGetBlockSize(coordinatesA, &d)); 1471 PetscCall(VecSetBlockSize(coordinatesB, d)); 1472 PetscCall(VecGetType(coordinatesA, &vtype)); 1473 PetscCall(VecSetType(coordinatesB, vtype)); 1474 PetscCall(VecGetArray(coordinatesA, &coordsA)); 1475 PetscCall(VecGetArray(coordinatesB, &coordsB)); 1476 for (v = 0; v < vEndB-vStartB; ++v) { 1477 PetscInt offA, offB; 1478 1479 PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA)); 1480 PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB)); 1481 for (d = 0; d < spaceDim; ++d) { 1482 coordsB[offB+d] = coordsA[offA+d]; 1483 } 1484 } 1485 if (lc) { /* localized coordinates */ 1486 PetscInt c; 1487 1488 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1489 PetscInt dof, offA, offB; 1490 1491 PetscCall(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA)); 1492 PetscCall(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB)); 1493 PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1494 PetscCall(PetscArraycpy(coordsB + offB,coordsA + offA,dof)); 1495 } 1496 } 1497 PetscCall(VecRestoreArray(coordinatesA, &coordsA)); 1498 PetscCall(VecRestoreArray(coordinatesB, &coordsB)); 1499 PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB)); 1500 PetscCall(VecDestroy(&coordinatesB)); 1501 PetscFunctionReturn(0); 1502 } 1503 1504 /*@ 1505 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1506 1507 Collective on dm 1508 1509 Input Parameter: 1510 . dm - The complete DMPlex object 1511 1512 Output Parameter: 1513 . dmUnint - The DMPlex object with only cells and vertices 1514 1515 Level: intermediate 1516 1517 Notes: 1518 It does not copy over the coordinates. 1519 1520 Developer Notes: 1521 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1522 1523 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1524 @*/ 1525 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1526 { 1527 DMPlexInterpolatedFlag interpolated; 1528 DM udm; 1529 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1530 1531 PetscFunctionBegin; 1532 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1533 PetscValidPointer(dmUnint, 2); 1534 PetscCall(DMGetDimension(dm, &dim)); 1535 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1536 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1537 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1538 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1539 PetscCall(PetscObjectReference((PetscObject) dm)); 1540 *dmUnint = dm; 1541 PetscFunctionReturn(0); 1542 } 1543 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1544 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1545 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), &udm)); 1546 PetscCall(DMSetType(udm, DMPLEX)); 1547 PetscCall(DMSetDimension(udm, dim)); 1548 PetscCall(DMPlexSetChart(udm, cStart, vEnd)); 1549 for (c = cStart; c < cEnd; ++c) { 1550 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1551 1552 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1553 for (cl = 0; cl < closureSize*2; cl += 2) { 1554 const PetscInt p = closure[cl]; 1555 1556 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1557 } 1558 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1559 PetscCall(DMPlexSetConeSize(udm, c, coneSize)); 1560 maxConeSize = PetscMax(maxConeSize, coneSize); 1561 } 1562 PetscCall(DMSetUp(udm)); 1563 PetscCall(PetscMalloc1(maxConeSize, &cone)); 1564 for (c = cStart; c < cEnd; ++c) { 1565 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1566 1567 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1568 for (cl = 0; cl < closureSize*2; cl += 2) { 1569 const PetscInt p = closure[cl]; 1570 1571 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1572 } 1573 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1574 PetscCall(DMPlexSetCone(udm, c, cone)); 1575 } 1576 PetscCall(PetscFree(cone)); 1577 PetscCall(DMPlexSymmetrize(udm)); 1578 PetscCall(DMPlexStratify(udm)); 1579 /* Reduce SF */ 1580 { 1581 PetscSF sfPoint, sfPointUn; 1582 const PetscSFNode *remotePoints; 1583 const PetscInt *localPoints; 1584 PetscSFNode *remotePointsUn; 1585 PetscInt *localPointsUn; 1586 PetscInt vEnd, numRoots, numLeaves, l; 1587 PetscInt numLeavesUn = 0, n = 0; 1588 1589 /* Get original SF information */ 1590 PetscCall(DMGetPointSF(dm, &sfPoint)); 1591 PetscCall(DMGetPointSF(udm, &sfPointUn)); 1592 PetscCall(DMPlexGetDepthStratum(dm, 0, NULL, &vEnd)); 1593 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1594 /* Allocate space for cells and vertices */ 1595 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1596 /* Fill in leaves */ 1597 if (vEnd >= 0) { 1598 PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn)); 1599 PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn)); 1600 for (l = 0; l < numLeaves; l++) { 1601 if (localPoints[l] < vEnd) { 1602 localPointsUn[n] = localPoints[l]; 1603 remotePointsUn[n].rank = remotePoints[l].rank; 1604 remotePointsUn[n].index = remotePoints[l].index; 1605 ++n; 1606 } 1607 } 1608 PetscCheck(n == numLeavesUn,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn); 1609 PetscCall(PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER)); 1610 } 1611 } 1612 /* This function makes the mesh fully uninterpolated on all ranks */ 1613 { 1614 DM_Plex *plex = (DM_Plex *) udm->data; 1615 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1616 } 1617 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, udm)); 1618 *dmUnint = udm; 1619 PetscFunctionReturn(0); 1620 } 1621 1622 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1623 { 1624 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1625 MPI_Comm comm; 1626 1627 PetscFunctionBegin; 1628 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1629 PetscCall(DMPlexGetDepth(dm, &depth)); 1630 PetscCall(DMGetDimension(dm, &dim)); 1631 1632 if (depth == dim) { 1633 *interpolated = DMPLEX_INTERPOLATED_FULL; 1634 if (!dim) goto finish; 1635 1636 /* Check points at height = dim are vertices (have no cones) */ 1637 PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd)); 1638 for (p=pStart; p<pEnd; p++) { 1639 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1640 if (coneSize) { 1641 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1642 goto finish; 1643 } 1644 } 1645 1646 /* Check points at height < dim have cones */ 1647 for (h=0; h<dim; h++) { 1648 PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 1649 for (p=pStart; p<pEnd; p++) { 1650 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1651 if (!coneSize) { 1652 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1653 goto finish; 1654 } 1655 } 1656 } 1657 } else if (depth == 1) { 1658 *interpolated = DMPLEX_INTERPOLATED_NONE; 1659 } else { 1660 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1661 } 1662 finish: 1663 PetscFunctionReturn(0); 1664 } 1665 1666 /*@ 1667 DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated. 1668 1669 Not Collective 1670 1671 Input Parameter: 1672 . dm - The DM object 1673 1674 Output Parameter: 1675 . interpolated - Flag whether the DM is interpolated 1676 1677 Level: intermediate 1678 1679 Notes: 1680 Unlike DMPlexIsInterpolatedCollective(), this is NOT collective 1681 so the results can be different on different ranks in special cases. 1682 However, DMPlexInterpolate() guarantees the result is the same on all. 1683 1684 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED. 1685 1686 Developer Notes: 1687 Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID. 1688 1689 If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called. 1690 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1691 DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL. 1692 1693 If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated. 1694 1695 DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL, 1696 and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1697 1698 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective() 1699 @*/ 1700 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1701 { 1702 DM_Plex *plex = (DM_Plex *) dm->data; 1703 1704 PetscFunctionBegin; 1705 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1706 PetscValidPointer(interpolated,2); 1707 if (plex->interpolated < 0) { 1708 PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated)); 1709 } else if (PetscDefined (USE_DEBUG)) { 1710 DMPlexInterpolatedFlag flg; 1711 1712 PetscCall(DMPlexIsInterpolated_Internal(dm, &flg)); 1713 PetscCheck(flg == plex->interpolated,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1714 } 1715 *interpolated = plex->interpolated; 1716 PetscFunctionReturn(0); 1717 } 1718 1719 /*@ 1720 DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner). 1721 1722 Collective 1723 1724 Input Parameter: 1725 . dm - The DM object 1726 1727 Output Parameter: 1728 . interpolated - Flag whether the DM is interpolated 1729 1730 Level: intermediate 1731 1732 Notes: 1733 Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks. 1734 1735 This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks. 1736 1737 Developer Notes: 1738 Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID. 1739 1740 If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated. 1741 MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1742 if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED, 1743 otherwise sets plex->interpolatedCollective = plex->interpolated. 1744 1745 If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective. 1746 1747 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated() 1748 @*/ 1749 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1750 { 1751 DM_Plex *plex = (DM_Plex *) dm->data; 1752 PetscBool debug=PETSC_FALSE; 1753 1754 PetscFunctionBegin; 1755 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1756 PetscValidPointer(interpolated,2); 1757 PetscCall(PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL)); 1758 if (plex->interpolatedCollective < 0) { 1759 DMPlexInterpolatedFlag min, max; 1760 MPI_Comm comm; 1761 1762 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1763 PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective)); 1764 PetscCallMPI(MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm)); 1765 PetscCallMPI(MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm)); 1766 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1767 if (debug) { 1768 PetscMPIInt rank; 1769 1770 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1771 PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective])); 1772 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 1773 } 1774 } 1775 *interpolated = plex->interpolatedCollective; 1776 PetscFunctionReturn(0); 1777 } 1778