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