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 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 66 if (faceTypes) CHKERRQ(DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &typesTmp)); 67 if (faceSizes) CHKERRQ(DMGetWorkArray(dm, PetscMax(maxConeSize, maxSupportSize), MPIU_INT, &sizesTmp)); 68 if (faces) CHKERRQ(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) CHKERRQ(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceTypes)); 321 if (faceSizes) CHKERRQ(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faceSizes)); 322 if (faces) CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 336 CHKERRQ(PetscHashIJKLCreate(&faceTable)); 337 CHKERRQ(PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES)); 338 CHKERRQ(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd)); 339 /* Number new faces and save face vertices in hash table */ 340 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, c, &ct)); 349 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 350 CHKERRQ(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 CHKERRQ(PetscSortInt(faceSize, (PetscInt *) &key)); 365 CHKERRQ(PetscHashIJKLPut(faceTable, key, &iter, &missing)); 366 if (missing) { 367 CHKERRQ(PetscHashIJKLIterSet(faceTable, iter, fEnd++)); 368 ++faceTypeNum[faceType]; 369 } 370 } 371 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, c, &ct)); 389 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 390 CHKERRQ(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 CHKERRQ(PetscSortInt(faceSize, (PetscInt *) &key)); 405 CHKERRQ(PetscHashIJKLPut(faceTable, key, &iter, &missing)); 406 if (missing) CHKERRQ(PetscHashIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++)); 407 } 408 CHKERRQ(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 409 } 410 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) { 411 PetscCheckFalse(faceTypeStart[ct] != faceTypeStart[ct-1] + faceTypeNum[ct],PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %D != %D + %D", DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct-1], faceTypeNum[ct]); 412 } 413 } 414 } 415 /* Add new points, always at the end of the numbering */ 416 CHKERRQ(DMPlexGetChart(dm, &pStart, &Np)); 417 CHKERRQ(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 CHKERRQ(DMCreateLabel(idm, "celltype")); 421 CHKERRQ(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 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 428 for (p = pStart; p < pEnd; ++p) { 429 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 430 CHKERRQ(DMPlexSetConeSize(idm, p, coneSize)); 431 CHKERRQ(DMPlexGetCellType(dm, p, &ct)); 432 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, c, &ct)); 442 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 443 CHKERRQ(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 444 CHKERRQ(DMPlexSetCellType(idm, c, ct)); 445 CHKERRQ(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 PetscCheckFalse(faceSize > 4,PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 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 CHKERRQ(PetscSortInt(faceSize, (PetscInt *) &key)); 461 CHKERRQ(PetscHashIJKLPut(faceTable, key, &iter, &missing)); 462 PetscCheckFalse(missing,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing face (cell %D, lf %D)", c, cf); 463 CHKERRQ(PetscHashIJKLIterGet(faceTable, iter, &f)); 464 CHKERRQ(DMPlexSetConeSize(idm, f, faceSize)); 465 CHKERRQ(DMPlexSetCellType(idm, f, faceType)); 466 } 467 CHKERRQ(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 468 } 469 CHKERRQ(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 CHKERRQ(DMPlexGetConeSection(idm, &cs)); 476 CHKERRQ(DMPlexGetCones(idm, &cones)); 477 CHKERRQ(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 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 487 for (p = pStart; p < pEnd; ++p) { 488 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 489 CHKERRQ(DMPlexSetCone(idm, p, cone)); 490 CHKERRQ(DMPlexGetConeOrientation(dm, p, &cone)); 491 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, c, &ct)); 501 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 502 CHKERRQ(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 PetscCheckFalse(faceSize > 4,PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %D > 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 CHKERRQ(PetscSortInt(faceSize, (PetscInt *) &key)); 519 CHKERRQ(PetscHashIJKLPut(faceTable, key, &iter, &missing)); 520 CHKERRQ(PetscHashIJKLIterGet(faceTable, iter, &f)); 521 CHKERRQ(DMPlexInsertCone(idm, c, cf, f)); 522 CHKERRQ(DMPlexGetCone(idm, f, &fcone)); 523 if (fcone[0] < 0) CHKERRQ(DMPlexSetCone(idm, f, face)); 524 { 525 const PetscInt *cone; 526 PetscInt coneSize, ornt; 527 528 CHKERRQ(DMPlexGetConeSize(idm, f, &coneSize)); 529 CHKERRQ(DMPlexGetCone(idm, f, &cone)); 530 PetscCheckFalse(coneSize != faceSize,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize); 531 /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */ 532 CHKERRQ(DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt)); 533 CHKERRQ(DMPlexInsertConeOrientation(idm, c, cf, ornt)); 534 } 535 } 536 CHKERRQ(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 537 } 538 CHKERRQ(PetscHashIJKLDestroy(&faceTable)); 539 CHKERRQ(DMPlexSymmetrize(idm)); 540 CHKERRQ(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 CHKERRQ(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote)); 554 nleaves = roffset[nranks]; 555 CHKERRQ(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 CHKERRQ(PetscArraycpy(&(*rmine1)[o], &rmine[o], n)); 562 CHKERRQ(PetscArraycpy(&(*rremote1)[o], &rremote[o], n)); 563 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 585 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 586 CHKERRMPI(MPI_Comm_size(comm, &size)); 587 CHKERRQ(DMGetPointSF(dm, &sf)); 588 CHKERRQ(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view")); 589 if (PetscDefined(USE_DEBUG)) CHKERRQ(DMPlexCheckPointSF(dm)); 590 CHKERRQ(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes)); 591 if (nroots < 0) PetscFunctionReturn(0); 592 CHKERRQ(PetscSFSetUp(sf)); 593 CHKERRQ(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1)); 594 for (p = 0; p < nleaves; ++p) { 595 PetscInt coneSize; 596 CHKERRQ(DMPlexGetConeSize(dm, locals[p], &coneSize)); 597 maxConeSize = PetscMax(maxConeSize, coneSize); 598 } 599 PetscCheckFalse(maxConeSize > 4,comm, PETSC_ERR_SUP, "This method does not support cones of size %D", maxConeSize); 600 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 606 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 639 CHKERRQ(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 640 CHKERRQ(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 641 CHKERRQ(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 642 if (debug) { 643 CHKERRQ(PetscSynchronizedFlush(comm, NULL)); 644 if (!rank) CHKERRQ(PetscSynchronizedPrintf(comm, "Referenced roots\n")); 645 } 646 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, p, &ct)); 654 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 655 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 656 if (debug) { 657 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] %4D: cone=[%4D %4D %4D %4D] roots=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D) (%d,%4D) (%d,%4D)]", 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 CHKERRQ(PetscFindMPIInt((PetscMPIInt) leavesRanks[p][c], nranks, ranks, &r)); 678 PetscCheckFalse(r < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leaf (%d,%D): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]); 679 PetscCheckFalse(ranks[r] < 0 || ranks[r] >= size,PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%D c=%D commsize=%d: ranks[%D] = %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 CHKERRQ(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0)); 684 PetscCheckFalse(ind0 < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): corresponding remote point not found - it seems there is missing connection in point SF!", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]); 685 /* Get the corresponding local point */ 686 mainCone[c] = rmine1[rS + ind0]; 687 } 688 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D %4D %4D]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3])); 689 /* Set the desired order of p's cone points and fix orientations accordingly */ 690 CHKERRQ(DMPolytopeGetOrientation(ct, cone, mainCone, &o)); 691 CHKERRQ(DMPlexOrientPoint(dm, p, o)); 692 } else if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, " ==\n")); 693 } 694 if (debug) { 695 CHKERRQ(PetscSynchronizedFlush(comm, NULL)); 696 CHKERRMPI(MPI_Barrier(comm)); 697 } 698 CHKERRQ(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view")); 699 CHKERRQ(PetscFree4(roots, leaves, rootsRanks, leavesRanks)); 700 CHKERRQ(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 CHKERRQ(PetscOptionsHasName(NULL, NULL, opt, &flg)); 712 if (!flg) PetscFunctionReturn(0); 713 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 714 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 715 for (idx = 0; idx < n; ++idx) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx])); 716 CHKERRQ(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 CHKERRQ(PetscOptionsHasName(NULL, NULL, opt, &flg)); 728 if (!flg) PetscFunctionReturn(0); 729 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 730 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 731 if (idxname) { 732 for (idx = 0; idx < n; ++idx) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d]%s %D rank %D index %D\n", rank, idxname, idx, a[idx].rank, a[idx].index)); 733 } else { 734 for (idx = 0; idx < n; ++idx) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index)); 735 } 736 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 748 CHKERRQ(DMGetPointSF(dm, &sf)); 749 CHKERRQ(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 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 778 CHKERRQ(DMGetPointSF(dm, &sf)); 779 CHKERRQ(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes)); 780 if (Nl < 0) goto owned; 781 CHKERRQ(PetscSFComputeDegreeBegin(sf, &rootdegree)); 782 CHKERRQ(PetscSFComputeDegreeEnd(sf, &rootdegree)); 783 if (rootdegree[localPoint]) goto owned; 784 CHKERRQ(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 CHKERRQ(DMGetPointSF(dm, &sf)); 803 CHKERRQ(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL)); 804 if (Nl < 0) PetscFunctionReturn(0); 805 CHKERRQ(PetscFindInt(p, Nl, locals, &idx)); 806 if (idx >= 0) {*isShared = PETSC_TRUE; PetscFunctionReturn(0);} 807 CHKERRQ(PetscSFComputeDegreeBegin(sf, &rootdegree)); 808 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 821 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 822 for (c = 0; c < coneSize; ++c) { 823 PetscBool pointShared; 824 825 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, p, &coneSize)); 840 CHKERRQ(DMPlexGetCone(dm, p, &cone)); 841 for (c = 0; c < coneSize; ++c) { 842 PetscSFNode rcp; 843 PetscBool mapFailed; 844 845 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 870 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 871 CHKERRQ(DMPlexGetOverlap(dm, &overlap)); 872 CHKERRQ(DMPlexGetVTKCellHeight(dm, &cellHeight)); 873 CHKERRQ(DMPlexGetPointHeight(dm, p, &height)); 874 if (!overlap && height <= cellHeight+1) { 875 /* cells can't be shared for non-overlapping meshes */ 876 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Skipping face %D to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p)); 877 PetscFunctionReturn(0); 878 } 879 CHKERRQ(DMPlexGetSupportSize(dm, p, &supportSize)); 880 CHKERRQ(DMPlexGetSupport(dm, p, &support)); 881 if (candidates) CHKERRQ(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) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face)); 892 key.i = p; 893 key.j = face; 894 CHKERRQ(PetscHMapIJGet(faceHash, key, &f)); 895 if (f >= 0) continue; 896 CHKERRQ(DMPlexConeIsShared(dm, face, &isShared)); 897 CHKERRQ(DMPlexGetConeMinimum(dm, face, &cpmin)); 898 CHKERRQ(DMPlexMapToGlobalPoint(dm, p, &rp, NULL)); 899 if (debug) { 900 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared)); 901 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Global point (%D, %D) Min Cone Point (%D, %D)\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index)); 902 } 903 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 904 CHKERRQ(PetscHMapIJSet(faceHash, key, p)); 905 if (candidates) { 906 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank)); 907 CHKERRQ(DMPlexGetConeSize(dm, face, &coneSize)); 908 CHKERRQ(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 CHKERRQ(DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx], NULL)); 918 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index)); 919 ++idx; 920 } 921 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "\n")); 922 } else { 923 /* Add cone size to section */ 924 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face)); 925 CHKERRQ(DMPlexGetConeSize(dm, face, &coneSize)); 926 CHKERRQ(PetscHMapIJSet(faceHash, key, p)); 927 CHKERRQ(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 CHKERRQ(DMPlexIsDistributed(dm, &flg)); 970 if (!flg) PetscFunctionReturn(0); 971 /* Set initial SF so that lower level queries work */ 972 CHKERRQ(DMSetPointSF(dm, pointSF)); 973 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 974 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 975 CHKERRQ(DMPlexGetOverlap(dm, &ov)); 976 PetscCheckFalse(ov,comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 977 CHKERRQ(PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug)); 978 CHKERRQ(PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view")); 979 CHKERRQ(PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view")); 980 CHKERRQ(PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0)); 981 /* Step 0: Precalculations */ 982 CHKERRQ(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints)); 983 PetscCheckFalse(Nr < 0,comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 984 CHKERRQ(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 CHKERRQ(PetscHMapIJSet(remoteHash, key, l)); 990 } 991 /* Compute root degree to identify shared points */ 992 CHKERRQ(PetscSFComputeDegreeBegin(pointSF, &rootdegree)); 993 CHKERRQ(PetscSFComputeDegreeEnd(pointSF, &rootdegree)); 994 CHKERRQ(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 CHKERRQ(PetscSectionCreate(comm, &candidateSection)); 1011 CHKERRQ(PetscSectionSetChart(candidateSection, 0, Nr)); 1012 { 1013 PetscHMapIJ faceHash; 1014 1015 CHKERRQ(PetscHMapIJCreate(&faceHash)); 1016 for (l = 0; l < Nl; ++l) { 1017 const PetscInt p = localPoints[l]; 1018 1019 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p)); 1020 CHKERRQ(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug)); 1021 } 1022 CHKERRQ(PetscHMapIJClear(faceHash)); 1023 CHKERRQ(PetscSectionSetUp(candidateSection)); 1024 CHKERRQ(PetscSectionGetStorageSize(candidateSection, &candidatesSize)); 1025 CHKERRQ(PetscMalloc1(candidatesSize, &candidates)); 1026 for (l = 0; l < Nl; ++l) { 1027 const PetscInt p = localPoints[l]; 1028 1029 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p)); 1030 CHKERRQ(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug)); 1031 } 1032 CHKERRQ(PetscHMapIJDestroy(&faceHash)); 1033 if (debug) CHKERRQ(PetscSynchronizedFlush(comm, NULL)); 1034 } 1035 CHKERRQ(PetscObjectSetName((PetscObject) candidateSection, "Candidate Section")); 1036 CHKERRQ(PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view")); 1037 CHKERRQ(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 CHKERRQ(PetscSFGetMultiSF(pointSF, &sfMulti)); 1045 CHKERRQ(PetscSFCreateInverseSF(sfMulti, &sfInverse)); 1046 CHKERRQ(PetscSectionCreate(comm, &candidateRemoteSection)); 1047 CHKERRQ(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection)); 1048 CHKERRQ(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates)); 1049 CHKERRQ(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize)); 1050 CHKERRQ(PetscMalloc1(candidatesRemoteSize, &candidatesRemote)); 1051 CHKERRQ(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE)); 1052 CHKERRQ(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE)); 1053 CHKERRQ(PetscSFDestroy(&sfInverse)); 1054 CHKERRQ(PetscSFDestroy(&sfCandidates)); 1055 CHKERRQ(PetscFree(remoteOffsets)); 1056 1057 CHKERRQ(PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section")); 1058 CHKERRQ(PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view")); 1059 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(candidateRemoteSection, idx, &dof)); 1075 CHKERRQ(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) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking face (%D, %D) at (%D, %D, %D) with cone size %D\n", rank, rface.rank, rface.index, r, idx, d, Np)); 1091 PetscCheckFalse(Np > 4,PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%D, %D) at (%D, %D, %D) with %D cone points", rface.rank, rface.index, r, idx, d, Np); 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 CHKERRQ(PetscSortSFNode(Np, (PetscSFNode *) &key)); 1101 CHKERRQ(PetscHashIJKLRemotePut(faceTable, key, &iter, &missing)); 1102 if (missing) { 1103 if (debug) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank)); 1104 CHKERRQ(PetscHashIJKLRemoteIterSet(faceTable, iter, rface)); 1105 } else { 1106 PetscSFNode oface; 1107 1108 CHKERRQ(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface)); 1109 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1110 if (debug) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank)); 1111 CHKERRQ(PetscHashIJKLRemoteIterSet(faceTable, iter, rface)); 1112 } 1113 } 1114 /* Check for local face */ 1115 points[0] = r; 1116 for (p = 1; p < Np; ++p) { 1117 CHKERRQ(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p], &mapToLocalPointFailed)); 1118 if (mapToLocalPointFailed) break; /* We got a point not in our overlap */ 1119 if (debug) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p])); 1120 } 1121 if (mapToLocalPointFailed) continue; 1122 CHKERRQ(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 CHKERRQ(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface)); 1131 if (debug) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing (%D, %D) with local face (%D, %D)\n", rank, oface.index, oface.rank, lface.index, lface.rank)); 1132 CHKERRQ(PetscHashIJKLRemoteIterSet(faceTable, iter, lface)); 1133 } 1134 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(candidateRemoteSection, idx2, &dof)); 1142 CHKERRQ(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) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx)); 1155 PetscCheckFalse(Np > 4,PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D 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 CHKERRQ(PetscSortSFNode(Np, (PetscSFNode *) &key)); 1165 if (debug) CHKERRQ(PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] key (%D, %D) (%D, %D) (%D, %D) (%D, %D)\n", rank, key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index)); 1166 CHKERRQ(PetscHashIJKLRemotePut(faceTable, key, &iter, &missing)); 1167 PetscCheckFalse(missing,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2); 1168 else CHKERRQ(PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx])); 1169 } 1170 } 1171 } 1172 if (debug) CHKERRQ(PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL)); 1173 CHKERRQ(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 CHKERRQ(PetscSFGetMultiSF(pointSF, &sfMulti)); 1184 CHKERRQ(PetscSectionCreate(comm, &claimSection)); 1185 CHKERRQ(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection)); 1186 CHKERRQ(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims)); 1187 CHKERRQ(PetscSectionGetStorageSize(claimSection, &claimsSize)); 1188 CHKERRQ(PetscMalloc1(claimsSize, &claims)); 1189 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1190 CHKERRQ(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE)); 1191 CHKERRQ(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE)); 1192 CHKERRQ(PetscSFDestroy(&sfClaims)); 1193 CHKERRQ(PetscFree(remoteOffsets)); 1194 CHKERRQ(PetscObjectSetName((PetscObject) claimSection, "Claim Section")); 1195 CHKERRQ(PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view")); 1196 CHKERRQ(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 CHKERRQ(PetscHMapICreate(&claimshash)); 1200 for (r = 0; r < Nr; ++r) { 1201 PetscInt dof, off, d; 1202 1203 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r)); 1204 CHKERRQ(PetscSectionGetDof(candidateSection, r, &dof)); 1205 CHKERRQ(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) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index)); 1214 points[0] = r; 1215 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0])); 1216 for (c = 0, d += 2; c < Np; ++c, ++d) { 1217 CHKERRQ(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1], NULL)); 1218 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1])); 1219 } 1220 CHKERRQ(DMPlexGetJoin(dm, Np+1, points, &joinSize, &join)); 1221 if (joinSize == 1) { 1222 if (claims[faceInd].rank == rank) { 1223 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0])); 1224 } else { 1225 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0])); 1226 CHKERRQ(PetscHMapISet(claimshash, join[0], faceInd)); 1227 } 1228 } else { 1229 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank)); 1230 } 1231 CHKERRQ(DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join)); 1232 } else { 1233 if (debug) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r)); 1234 d += claims[off+d].index+1; 1235 } 1236 } 1237 } 1238 if (debug) CHKERRQ(PetscSynchronizedFlush(comm, NULL)); 1239 /* Step 6) Create new pointSF from hashed claims */ 1240 CHKERRQ(PetscHMapIGetSize(claimshash, &NlNew)); 1241 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 1242 CHKERRQ(PetscMalloc1(Nl + NlNew, &localPointsNew)); 1243 CHKERRQ(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 CHKERRQ(PetscHMapIGetKeys(claimshash, &p, localPointsNew)); 1251 /* We sort new points, and assume they are numbered after all existing points */ 1252 CHKERRQ(PetscSortInt(NlNew, &localPointsNew[Nl])); 1253 for (p = Nl; p < Nl + NlNew; ++p) { 1254 PetscInt off; 1255 CHKERRQ(PetscHMapIGet(claimshash, localPointsNew[p], &off)); 1256 PetscCheckFalse(claims[off].rank < 0 || claims[off].index < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %D, (%D, %D)", localPointsNew[p], claims[off].rank, claims[off].index); 1257 remotePointsNew[p] = claims[off]; 1258 } 1259 CHKERRQ(PetscSFCreate(comm, &sfPointNew)); 1260 CHKERRQ(PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1261 CHKERRQ(PetscSFSetUp(sfPointNew)); 1262 CHKERRQ(DMSetPointSF(dm, sfPointNew)); 1263 CHKERRQ(PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view")); 1264 CHKERRQ(PetscSFDestroy(&sfPointNew)); 1265 CHKERRQ(PetscHMapIDestroy(&claimshash)); 1266 } 1267 CHKERRQ(PetscHMapIJDestroy(&remoteHash)); 1268 CHKERRQ(PetscSectionDestroy(&candidateSection)); 1269 CHKERRQ(PetscSectionDestroy(&candidateRemoteSection)); 1270 CHKERRQ(PetscSectionDestroy(&claimSection)); 1271 CHKERRQ(PetscFree(candidates)); 1272 CHKERRQ(PetscFree(candidatesRemote)); 1273 CHKERRQ(PetscFree(claims)); 1274 CHKERRQ(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 It does not copy over the coordinates. 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 CHKERRQ(PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0)); 1313 CHKERRQ(DMPlexGetDepth(dm, &depth)); 1314 CHKERRQ(DMGetDimension(dm, &dim)); 1315 CHKERRQ(DMPlexIsInterpolated(dm, &interpolated)); 1316 PetscCheckFalse(interpolated == DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1317 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1318 CHKERRQ(PetscObjectReference((PetscObject) dm)); 1319 idm = dm; 1320 } else { 1321 for (d = 1; d < dim; ++d) { 1322 /* Create interpolated mesh */ 1323 CHKERRQ(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1324 CHKERRQ(DMSetType(idm, DMPLEX)); 1325 CHKERRQ(DMSetDimension(idm, dim)); 1326 if (depth > 0) { 1327 CHKERRQ(DMPlexInterpolateFaces_Internal(odm, 1, idm)); 1328 CHKERRQ(DMGetPointSF(odm, &sfPoint)); 1329 { 1330 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1331 PetscInt nroots; 1332 CHKERRQ(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1333 if (nroots >= 0) CHKERRQ(DMPlexInterpolatePointSF(idm, sfPoint)); 1334 } 1335 } 1336 if (odm != dm) CHKERRQ(DMDestroy(&odm)); 1337 odm = idm; 1338 } 1339 CHKERRQ(PetscObjectGetName((PetscObject) dm, &name)); 1340 CHKERRQ(PetscObjectSetName((PetscObject) idm, name)); 1341 CHKERRQ(DMPlexCopyCoordinates(dm, idm)); 1342 CHKERRQ(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1343 CHKERRQ(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL)); 1344 if (flg) CHKERRQ(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 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, idm)); 1352 *dmInt = idm; 1353 CHKERRQ(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 CHKERRQ(DMGetCoordinateDim(dmA, &cdim)); 1389 CHKERRQ(DMSetCoordinateDim(dmB, cdim)); 1390 CHKERRQ(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA)); 1391 CHKERRQ(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB)); 1392 PetscCheckFalse((vEndA-vStartA) != (vEndB-vStartB),PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d 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 CHKERRQ(DMGetCoordinateDM(dmA, &cdmA)); 1403 CHKERRQ(DMGetCoordinateDM(dmB, &cdmB)); 1404 CHKERRQ(DMGetField(cdmA, 0, NULL, &objA)); 1405 CHKERRQ(DMGetField(cdmB, 0, NULL, &objB)); 1406 CHKERRQ(PetscObjectGetClassId(objA, &idA)); 1407 CHKERRQ(PetscObjectGetClassId(objB, &idB)); 1408 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1409 CHKERRQ(DMSetField(cdmB, 0, NULL, objA)); 1410 CHKERRQ(DMCreateDS(cdmB)); 1411 CHKERRQ(DMGetDS(cdmA, &dsA)); 1412 CHKERRQ(DMGetDS(cdmB, &dsB)); 1413 CHKERRQ(PetscDSGetCoordinateDimension(dsA, &cdim)); 1414 CHKERRQ(PetscDSSetCoordinateDimension(dsB, cdim)); 1415 CHKERRQ(PetscDSGetConstants(dsA, &Nc, &constants)); 1416 CHKERRQ(PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants)); 1417 } 1418 } 1419 CHKERRQ(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA)); 1420 CHKERRQ(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB)); 1421 CHKERRQ(DMGetCoordinateSection(dmA, &coordSectionA)); 1422 CHKERRQ(DMGetCoordinateSection(dmB, &coordSectionB)); 1423 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1424 CHKERRQ(PetscSectionGetNumFields(coordSectionA, &Nf)); 1425 if (!Nf) PetscFunctionReturn(0); 1426 PetscCheckFalse(Nf > 1,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1427 if (!coordSectionB) { 1428 PetscInt dim; 1429 1430 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB)); 1431 CHKERRQ(DMGetCoordinateDim(dmA, &dim)); 1432 CHKERRQ(DMSetCoordinateSection(dmB, dim, coordSectionB)); 1433 CHKERRQ(PetscObjectDereference((PetscObject) coordSectionB)); 1434 } 1435 CHKERRQ(PetscSectionSetNumFields(coordSectionB, 1)); 1436 CHKERRQ(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim)); 1437 CHKERRQ(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim)); 1438 CHKERRQ(PetscSectionGetChart(coordSectionA, &cS, &cE)); 1439 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1440 PetscCheckFalse((cEndA-cStartA) != (cEndB-cStartB),PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D 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 CHKERRQ(PetscSectionSetChart(coordSectionB, cS, cE)); 1449 for (v = vStartB; v < vEndB; ++v) { 1450 CHKERRQ(PetscSectionSetDof(coordSectionB, v, spaceDim)); 1451 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1460 CHKERRQ(PetscSectionSetDof(coordSectionB, c + cStartB, dof)); 1461 CHKERRQ(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof)); 1462 } 1463 } 1464 CHKERRQ(PetscSectionSetUp(coordSectionB)); 1465 CHKERRQ(PetscSectionGetStorageSize(coordSectionB, &coordSizeB)); 1466 CHKERRQ(DMGetCoordinatesLocal(dmA, &coordinatesA)); 1467 CHKERRQ(VecCreate(PETSC_COMM_SELF, &coordinatesB)); 1468 CHKERRQ(PetscObjectSetName((PetscObject) coordinatesB, "coordinates")); 1469 CHKERRQ(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE)); 1470 CHKERRQ(VecGetBlockSize(coordinatesA, &d)); 1471 CHKERRQ(VecSetBlockSize(coordinatesB, d)); 1472 CHKERRQ(VecGetType(coordinatesA, &vtype)); 1473 CHKERRQ(VecSetType(coordinatesB, vtype)); 1474 CHKERRQ(VecGetArray(coordinatesA, &coordsA)); 1475 CHKERRQ(VecGetArray(coordinatesB, &coordsB)); 1476 for (v = 0; v < vEndB-vStartB; ++v) { 1477 PetscInt offA, offB; 1478 1479 CHKERRQ(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA)); 1480 CHKERRQ(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 CHKERRQ(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA)); 1492 CHKERRQ(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB)); 1493 CHKERRQ(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1494 CHKERRQ(PetscArraycpy(coordsB + offB,coordsA + offA,dof)); 1495 } 1496 } 1497 CHKERRQ(VecRestoreArray(coordinatesA, &coordsA)); 1498 CHKERRQ(VecRestoreArray(coordinatesB, &coordsB)); 1499 CHKERRQ(DMSetCoordinatesLocal(dmB, coordinatesB)); 1500 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 1535 CHKERRQ(DMPlexIsInterpolated(dm, &interpolated)); 1536 PetscCheckFalse(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 CHKERRQ(PetscObjectReference((PetscObject) dm)); 1540 *dmUnint = dm; 1541 PetscFunctionReturn(0); 1542 } 1543 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1544 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1545 CHKERRQ(DMCreate(PetscObjectComm((PetscObject) dm), &udm)); 1546 CHKERRQ(DMSetType(udm, DMPLEX)); 1547 CHKERRQ(DMSetDimension(udm, dim)); 1548 CHKERRQ(DMPlexSetChart(udm, cStart, vEnd)); 1549 for (c = cStart; c < cEnd; ++c) { 1550 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1551 1552 CHKERRQ(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 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1559 CHKERRQ(DMPlexSetConeSize(udm, c, coneSize)); 1560 maxConeSize = PetscMax(maxConeSize, coneSize); 1561 } 1562 CHKERRQ(DMSetUp(udm)); 1563 CHKERRQ(PetscMalloc1(maxConeSize, &cone)); 1564 for (c = cStart; c < cEnd; ++c) { 1565 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1566 1567 CHKERRQ(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 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1574 CHKERRQ(DMPlexSetCone(udm, c, cone)); 1575 } 1576 CHKERRQ(PetscFree(cone)); 1577 CHKERRQ(DMPlexSymmetrize(udm)); 1578 CHKERRQ(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 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 1591 CHKERRQ(DMGetPointSF(udm, &sfPointUn)); 1592 CHKERRQ(DMPlexGetDepthStratum(dm, 0, NULL, &vEnd)); 1593 CHKERRQ(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 CHKERRQ(PetscMalloc1(numLeavesUn, &remotePointsUn)); 1599 CHKERRQ(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 PetscCheckFalse(n != numLeavesUn,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1609 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1629 CHKERRQ(DMPlexGetDepth(dm, &depth)); 1630 CHKERRQ(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 CHKERRQ(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd)); 1638 for (p=pStart; p<pEnd; p++) { 1639 CHKERRQ(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 CHKERRQ(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 1649 for (p=pStart; p<pEnd; p++) { 1650 CHKERRQ(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 CHKERRQ(DMPlexIsInterpolated_Internal(dm, &plex->interpolated)); 1709 } else if (PetscDefined (USE_DEBUG)) { 1710 DMPlexInterpolatedFlag flg; 1711 1712 CHKERRQ(DMPlexIsInterpolated_Internal(dm, &flg)); 1713 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 1763 CHKERRQ(DMPlexIsInterpolated(dm, &plex->interpolatedCollective)); 1764 CHKERRMPI(MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm)); 1765 CHKERRMPI(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 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1771 CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective])); 1772 CHKERRQ(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 1773 } 1774 } 1775 *interpolated = plex->interpolatedCollective; 1776 PetscFunctionReturn(0); 1777 } 1778