1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/hashmapi.h> 3 #include <petsc/private/hashmapij.h> 4 5 const char * const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL}; 6 7 /* HashIJKL */ 8 9 #include <petsc/private/hashmap.h> 10 11 typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey; 12 13 #define PetscHashIJKLKeyHash(key) \ 14 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \ 15 PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l))) 16 17 #define PetscHashIJKLKeyEqual(k1,k2) \ 18 (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0) 19 20 PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1) 21 22 static PetscSFNode _PetscInvalidSFNode = {-1, -1}; 23 24 typedef struct _PetscHashIJKLRemoteKey { PetscSFNode i, j, k, l; } PetscHashIJKLRemoteKey; 25 26 #define PetscHashIJKLRemoteKeyHash(key) \ 27 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index),PetscHashInt((key).j.rank + (key).j.index)), \ 28 PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index),PetscHashInt((key).l.rank + (key).l.index))) 29 30 #define PetscHashIJKLRemoteKeyEqual(k1,k2) \ 31 (((k1).i.rank==(k2).i.rank) ? ((k1).i.index==(k2).i.index) ? ((k1).j.rank==(k2).j.rank) ? ((k1).j.index==(k2).j.index) ? ((k1).k.rank==(k2).k.rank) ? ((k1).k.index==(k2).k.index) ? ((k1).l.rank==(k2).l.rank) ? ((k1).l.index==(k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0) 32 33 PETSC_HASH_MAP(HashIJKLRemote, PetscHashIJKLRemoteKey, PetscSFNode, PetscHashIJKLRemoteKeyHash, PetscHashIJKLRemoteKeyEqual, _PetscInvalidSFNode) 34 35 static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[]) 36 { 37 PetscInt i; 38 39 PetscFunctionBegin; 40 for (i = 1; i < n; ++i) { 41 PetscSFNode x = A[i]; 42 PetscInt j; 43 44 for (j = i-1; j >= 0; --j) { 45 if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break; 46 A[j+1] = A[j]; 47 } 48 A[j+1] = x; 49 } 50 PetscFunctionReturn(0); 51 } 52 53 /* 54 DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone 55 */ 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: SETERRQ1(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, 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 if (faceSize > 4) SETERRQ1(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 if (faceSize > 4) SETERRQ1(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 if (faceTypeStart[ct] != faceTypeStart[ct-1] + faceTypeNum[ct]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %D != %D + %D", DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct-1], faceTypeNum[ct]); 416 } 417 } 418 } 419 /* Add new points, always at the end of the numbering */ 420 ierr = DMPlexGetChart(dm, NULL, &Np);CHKERRQ(ierr); 421 ierr = DMPlexSetChart(idm, 0, 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 if (faceSize > 4) SETERRQ1(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 if (missing) SETERRQ2(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 if (faceSize > 4) SETERRQ1(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 /* TODO This should be unnecessary since we have autoamtic orientation */ 529 { 530 /* when matching hybrid faces in 3D, only few cases are possible. 531 Face traversal however can no longer follow the usual convention, this seems a serious issue to me */ 532 PetscInt tquad_map[4][4] = { {PETSC_MIN_INT, 0,PETSC_MIN_INT,PETSC_MIN_INT}, 533 { -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT}, 534 {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT, 1}, 535 {PETSC_MIN_INT,PETSC_MIN_INT, -2,PETSC_MIN_INT} }; 536 PetscInt i, i2, j; 537 const PetscInt *cone; 538 PetscInt coneSize, ornt; 539 540 /* Orient face: Do not allow reverse orientation at the first vertex */ 541 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 542 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 543 if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize); 544 /* - First find the initial vertex */ 545 for (i = 0; i < faceSize; ++i) if (face[0] == cone[i]) break; 546 /* If we want to compare tensor faces to regular faces, we have to flip them and take the else branch here */ 547 if (faceType == DM_POLYTOPE_SEG_PRISM_TENSOR) { 548 /* find the second vertex */ 549 for (i2 = 0; i2 < faceSize; ++i2) if (face[1] == cone[i2]) break; 550 switch (faceSize) { 551 case 2: 552 ornt = i ? -2 : 0; 553 break; 554 case 4: 555 ornt = tquad_map[i][i2]; 556 break; 557 default: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSize, f, c); 558 } 559 } else { 560 /* Try forward comparison */ 561 for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+j)%faceSize]) break; 562 if (j == faceSize) { 563 if ((faceSize == 2) && (i == 1)) ornt = -2; 564 else ornt = i; 565 } else { 566 /* Try backward comparison */ 567 for (j = 0; j < faceSize; ++j) if (face[j] != cone[(i+faceSize-j)%faceSize]) break; 568 if (j == faceSize) { 569 if (i == 0) ornt = -faceSize; 570 else ornt = -i; 571 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c); 572 } 573 } 574 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 575 } 576 } 577 ierr = DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces);CHKERRQ(ierr); 578 } 579 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 580 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 581 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 582 PetscFunctionReturn(0); 583 } 584 585 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 586 { 587 PetscInt nleaves; 588 PetscInt nranks; 589 const PetscMPIInt *ranks=NULL; 590 const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL; 591 PetscInt n, o, r; 592 PetscErrorCode ierr; 593 594 PetscFunctionBegin; 595 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 596 nleaves = roffset[nranks]; 597 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr); 598 for (r=0; r<nranks; r++) { 599 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 600 - to unify order with the other side */ 601 o = roffset[r]; 602 n = roffset[r+1] - o; 603 ierr = PetscArraycpy(&(*rmine1)[o], &rmine[o], n);CHKERRQ(ierr); 604 ierr = PetscArraycpy(&(*rremote1)[o], &rremote[o], n);CHKERRQ(ierr); 605 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr); 606 } 607 PetscFunctionReturn(0); 608 } 609 610 PetscErrorCode DMPlexOrientInterface_Internal(DM dm) 611 { 612 /* Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */ 613 PetscInt masterCone[2]; 614 PetscInt (*roots)[2], (*leaves)[2]; 615 PetscMPIInt (*rootsRanks)[2], (*leavesRanks)[2]; 616 617 PetscSF sf=NULL; 618 const PetscInt *locals=NULL; 619 const PetscSFNode *remotes=NULL; 620 PetscInt nroots, nleaves, p, c; 621 PetscInt nranks, n, o, r; 622 const PetscMPIInt *ranks=NULL; 623 const PetscInt *roffset=NULL; 624 PetscInt *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 625 const PetscInt *cone=NULL; 626 PetscInt coneSize, ind0; 627 MPI_Comm comm; 628 PetscMPIInt rank,size; 629 PetscInt debug = 0; 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 634 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);CHKERRQ(ierr); 635 if (nroots < 0) PetscFunctionReturn(0); 636 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 637 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr); 638 ierr = DMViewFromOptions(dm, NULL, "-before_fix_dm_view");CHKERRQ(ierr); 639 if (PetscDefined(USE_DEBUG)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 640 ierr = SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr); 641 ierr = PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);CHKERRQ(ierr); 642 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 643 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 644 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 645 for (p = 0; p < nroots; ++p) { 646 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 647 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 648 if (coneSize < 2) { 649 for (c = 0; c < 2; c++) { 650 roots[p][c] = -1; 651 rootsRanks[p][c] = -1; 652 } 653 continue; 654 } 655 /* Translate all points to root numbering */ 656 for (c = 0; c < 2; c++) { 657 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr); 658 if (ind0 < 0) { 659 roots[p][c] = cone[c]; 660 rootsRanks[p][c] = rank; 661 } else { 662 roots[p][c] = remotes[ind0].index; 663 rootsRanks[p][c] = remotes[ind0].rank; 664 } 665 } 666 } 667 for (p = 0; p < nroots; ++p) { 668 for (c = 0; c < 2; c++) { 669 leaves[p][c] = -2; 670 leavesRanks[p][c] = -2; 671 } 672 } 673 ierr = PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 674 ierr = PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 675 ierr = PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 676 ierr = PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 677 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 678 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Referenced roots\n");CHKERRQ(ierr);} 679 for (p = 0; p < nroots; ++p) { 680 if (leaves[p][0] < 0) continue; 681 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 682 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 683 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] %4D: cone=[%4D %4D] roots=[(%d,%4D) (%d,%4D)] leaves=[(%d,%4D) (%d,%4D)]", rank, p, cone[0], cone[1], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1]);CHKERRQ(ierr);} 684 if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][1] != rootsRanks[p][1])) { 685 /* Translate these two leaves to my cone points; masterCone means desired order p's cone points */ 686 for (c = 0; c < 2; c++) { 687 if (leavesRanks[p][c] == rank) { 688 /* A local leave is just taken as it is */ 689 masterCone[c] = leaves[p][c]; 690 continue; 691 } 692 /* Find index of rank leavesRanks[p][c] among remote ranks */ 693 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */ 694 ierr = PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);CHKERRQ(ierr); 695 if (PetscUnlikely(r < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): leave rank not found among remote ranks",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]); 696 if (PetscUnlikely(ranks[r] < 0 || ranks[r] >= size)) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%D c=%D commsize=%d: ranks[%D] = %d makes no sense",p,c,size,r,ranks[r]); 697 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */ 698 o = roffset[r]; 699 n = roffset[r+1] - o; 700 ierr = PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);CHKERRQ(ierr); 701 if (PetscUnlikely(ind0 < 0)) SETERRQ7(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D cone[%D]=%D root (%d,%D) leave (%d,%D): corresponding remote point not found - it seems there is missing connection in point SF!",p,c,cone[c],rootsRanks[p][c],roots[p][c],leavesRanks[p][c],leaves[p][c]); 702 /* Get the corresponding local point */ 703 masterCone[c] = rmine1[o+ind0];CHKERRQ(ierr); 704 } 705 if (debug) {ierr = PetscSynchronizedPrintf(comm, " masterCone=[%4D %4D]\n", masterCone[0], masterCone[1]);CHKERRQ(ierr);} 706 /* Set the desired order of p's cone points and fix orientations accordingly */ 707 /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */ 708 ierr = DMPlexOrientCell(dm, p, 2, masterCone);CHKERRQ(ierr); 709 } else if (debug) {ierr = PetscSynchronizedPrintf(comm, " ==\n");CHKERRQ(ierr);} 710 } 711 if (debug) { 712 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 713 ierr = MPI_Barrier(comm);CHKERRQ(ierr); 714 } 715 ierr = DMViewFromOptions(dm, NULL, "-after_fix_dm_view");CHKERRQ(ierr); 716 ierr = PetscFree4(roots, leaves, rootsRanks, leavesRanks);CHKERRQ(ierr); 717 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr); 718 PetscFunctionReturn(0); 719 } 720 721 static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[]) 722 { 723 PetscInt idx; 724 PetscMPIInt rank; 725 PetscBool flg; 726 PetscErrorCode ierr; 727 728 PetscFunctionBegin; 729 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr); 730 if (!flg) PetscFunctionReturn(0); 731 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 732 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr); 733 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);CHKERRQ(ierr);} 734 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 735 PetscFunctionReturn(0); 736 } 737 738 static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[]) 739 { 740 PetscInt idx; 741 PetscMPIInt rank; 742 PetscBool flg; 743 PetscErrorCode ierr; 744 745 PetscFunctionBegin; 746 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr); 747 if (!flg) PetscFunctionReturn(0); 748 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 749 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr); 750 if (idxname) { 751 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);} 752 } else { 753 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);CHKERRQ(ierr);} 754 } 755 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 756 PetscFunctionReturn(0); 757 } 758 759 static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint) 760 { 761 PetscSF sf; 762 const PetscInt *locals; 763 PetscMPIInt rank; 764 PetscErrorCode ierr; 765 766 PetscFunctionBegin; 767 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 768 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 769 ierr = PetscSFGetGraph(sf, NULL, NULL, &locals, NULL);CHKERRQ(ierr); 770 if (remotePoint.rank == rank) { 771 *localPoint = remotePoint.index; 772 } else { 773 PetscHashIJKey key; 774 PetscInt l; 775 776 key.i = remotePoint.index; 777 key.j = remotePoint.rank; 778 ierr = PetscHMapIJGet(remotehash, key, &l);CHKERRQ(ierr); 779 if (l >= 0) { 780 *localPoint = locals[l]; 781 } else PetscFunctionReturn(1); 782 } 783 PetscFunctionReturn(0); 784 } 785 786 static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint) 787 { 788 PetscSF sf; 789 const PetscInt *locals, *rootdegree; 790 const PetscSFNode *remotes; 791 PetscInt Nl, l; 792 PetscMPIInt rank; 793 PetscErrorCode ierr; 794 795 PetscFunctionBegin; 796 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 797 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 798 ierr = PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes);CHKERRQ(ierr); 799 if (Nl < 0) goto owned; 800 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 801 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 802 if (rootdegree[localPoint]) goto owned; 803 ierr = PetscFindInt(localPoint, Nl, locals, &l);CHKERRQ(ierr); 804 if (l < 0) PetscFunctionReturn(1); 805 *remotePoint = remotes[l]; 806 PetscFunctionReturn(0); 807 owned: 808 remotePoint->rank = rank; 809 remotePoint->index = localPoint; 810 PetscFunctionReturn(0); 811 } 812 813 814 static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared) 815 { 816 PetscSF sf; 817 const PetscInt *locals, *rootdegree; 818 PetscInt Nl, idx; 819 PetscErrorCode ierr; 820 821 PetscFunctionBegin; 822 *isShared = PETSC_FALSE; 823 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 824 ierr = PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);CHKERRQ(ierr); 825 if (Nl < 0) PetscFunctionReturn(0); 826 ierr = PetscFindInt(p, Nl, locals, &idx);CHKERRQ(ierr); 827 if (idx >= 0) {*isShared = PETSC_TRUE; PetscFunctionReturn(0);} 828 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 829 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 830 if (rootdegree[p] > 0) *isShared = PETSC_TRUE; 831 PetscFunctionReturn(0); 832 } 833 834 static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared) 835 { 836 const PetscInt *cone; 837 PetscInt coneSize, c; 838 PetscBool cShared = PETSC_TRUE; 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 843 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 844 for (c = 0; c < coneSize; ++c) { 845 PetscBool pointShared; 846 847 ierr = DMPlexPointIsShared(dm, cone[c], &pointShared);CHKERRQ(ierr); 848 cShared = (PetscBool) (cShared && pointShared); 849 } 850 *isShared = coneSize ? cShared : PETSC_FALSE; 851 PetscFunctionReturn(0); 852 } 853 854 static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin) 855 { 856 const PetscInt *cone; 857 PetscInt coneSize, c; 858 PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1}; 859 PetscErrorCode ierr; 860 861 PetscFunctionBegin; 862 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 863 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 864 for (c = 0; c < coneSize; ++c) { 865 PetscSFNode rcp; 866 867 ierr = DMPlexMapToGlobalPoint(dm, cone[c], &rcp); 868 if (ierr) { 869 cmin = missing; 870 } else { 871 cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin; 872 } 873 } 874 *cpmin = coneSize ? cmin : missing; 875 PetscFunctionReturn(0); 876 } 877 878 /* 879 Each shared face has an entry in the candidates array: 880 (-1, coneSize-1), {(global cone point)} 881 where the set is missing the point p which we use as the key for the face 882 */ 883 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug) 884 { 885 MPI_Comm comm; 886 const PetscInt *support; 887 PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height; 888 PetscMPIInt rank; 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 893 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 894 ierr = DMPlexGetOverlap(dm, &overlap);CHKERRQ(ierr); 895 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 896 ierr = DMPlexGetPointHeight(dm, p, &height);CHKERRQ(ierr); 897 if (!overlap && height <= cellHeight+1) { 898 /* cells can't be shared for non-overlapping meshes */ 899 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);} 900 PetscFunctionReturn(0); 901 } 902 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 903 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 904 if (candidates) {ierr = PetscSectionGetOffset(candidateSection, p, &off);CHKERRQ(ierr);} 905 for (s = 0; s < supportSize; ++s) { 906 const PetscInt face = support[s]; 907 const PetscInt *cone; 908 PetscSFNode cpmin={-1,-1}, rp={-1,-1}; 909 PetscInt coneSize, c, f; 910 PetscBool isShared = PETSC_FALSE; 911 PetscHashIJKey key; 912 913 /* Only add point once */ 914 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face);CHKERRQ(ierr);} 915 key.i = p; 916 key.j = face; 917 ierr = PetscHMapIJGet(faceHash, key, &f);CHKERRQ(ierr); 918 if (f >= 0) continue; 919 ierr = DMPlexConeIsShared(dm, face, &isShared);CHKERRQ(ierr); 920 ierr = DMPlexGetConeMinimum(dm, face, &cpmin);CHKERRQ(ierr); 921 ierr = DMPlexMapToGlobalPoint(dm, p, &rp);CHKERRQ(ierr); 922 if (debug) { 923 ierr = PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared);CHKERRQ(ierr); 924 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); 925 } 926 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 927 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 928 if (candidates) { 929 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank);CHKERRQ(ierr);} 930 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 931 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 932 candidates[off+idx].rank = -1; 933 candidates[off+idx++].index = coneSize-1; 934 candidates[off+idx].rank = rank; 935 candidates[off+idx++].index = face; 936 for (c = 0; c < coneSize; ++c) { 937 const PetscInt cp = cone[c]; 938 939 if (cp == p) continue; 940 ierr = DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);CHKERRQ(ierr); 941 if (debug) {ierr = PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);CHKERRQ(ierr);} 942 ++idx; 943 } 944 if (debug) {ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr);} 945 } else { 946 /* Add cone size to section */ 947 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face);CHKERRQ(ierr);} 948 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 949 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 950 ierr = PetscSectionAddDof(candidateSection, p, coneSize+1);CHKERRQ(ierr); 951 } 952 } 953 } 954 PetscFunctionReturn(0); 955 } 956 957 /*@ 958 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation 959 960 Collective on dm 961 962 Input Parameters: 963 + dm - The interpolated DM 964 - pointSF - The initial SF without interpolated points 965 966 Output Parameter: 967 . pointSF - The SF including interpolated points 968 969 Level: developer 970 971 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 972 973 .seealso: DMPlexInterpolate(), DMPlexUninterpolate() 974 @*/ 975 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 976 { 977 MPI_Comm comm; 978 PetscHMapIJ remoteHash; 979 PetscHMapI claimshash; 980 PetscSection candidateSection, candidateRemoteSection, claimSection; 981 PetscSFNode *candidates, *candidatesRemote, *claims; 982 const PetscInt *localPoints, *rootdegree; 983 const PetscSFNode *remotePoints; 984 PetscInt ov, Nr, r, Nl, l; 985 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 986 PetscBool flg, debug = PETSC_FALSE; 987 PetscMPIInt rank; 988 PetscErrorCode ierr; 989 990 PetscFunctionBegin; 991 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 992 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 3); 993 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 994 if (!flg) PetscFunctionReturn(0); 995 /* Set initial SF so that lower level queries work */ 996 ierr = DMSetPointSF(dm, pointSF);CHKERRQ(ierr); 997 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 998 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 999 ierr = DMPlexGetOverlap(dm, &ov);CHKERRQ(ierr); 1000 if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1001 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr); 1002 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr); 1003 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr); 1004 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1005 /* Step 0: Precalculations */ 1006 ierr = PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);CHKERRQ(ierr); 1007 if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1008 ierr = PetscHMapIJCreate(&remoteHash);CHKERRQ(ierr); 1009 for (l = 0; l < Nl; ++l) { 1010 PetscHashIJKey key; 1011 key.i = remotePoints[l].index; 1012 key.j = remotePoints[l].rank; 1013 ierr = PetscHMapIJSet(remoteHash, key, l);CHKERRQ(ierr); 1014 } 1015 /* Compute root degree to identify shared points */ 1016 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 1017 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 1018 ierr = IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);CHKERRQ(ierr); 1019 /* 1020 1) Loop over each leaf point $p$ at depth $d$ in the SF 1021 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1022 \begin{itemize} 1023 \item all cone points of $f$ are shared 1024 \item $p$ is the cone point with smallest canonical number 1025 \end{itemize} 1026 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1027 \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 1028 \item Send the root face from the root back to all leaf process 1029 \item Leaf processes add the shared face to the SF 1030 */ 1031 /* Step 1: Construct section+SFNode array 1032 The section has entries for all shared faces for which we have a leaf point in the cone 1033 The array holds candidate shared faces, each face is refered to by the leaf point */ 1034 ierr = PetscSectionCreate(comm, &candidateSection);CHKERRQ(ierr); 1035 ierr = PetscSectionSetChart(candidateSection, 0, Nr);CHKERRQ(ierr); 1036 { 1037 PetscHMapIJ faceHash; 1038 1039 ierr = PetscHMapIJCreate(&faceHash);CHKERRQ(ierr); 1040 for (l = 0; l < Nl; ++l) { 1041 const PetscInt p = localPoints[l]; 1042 1043 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1044 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);CHKERRQ(ierr); 1045 } 1046 ierr = PetscHMapIJClear(faceHash);CHKERRQ(ierr); 1047 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 1048 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 1049 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 1050 for (l = 0; l < Nl; ++l) { 1051 const PetscInt p = localPoints[l]; 1052 1053 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1054 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);CHKERRQ(ierr); 1055 } 1056 ierr = PetscHMapIJDestroy(&faceHash);CHKERRQ(ierr); 1057 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1058 } 1059 ierr = PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");CHKERRQ(ierr); 1060 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr); 1061 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr); 1062 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1063 /* Note that this section is indexed by offsets into leaves, not by point number */ 1064 { 1065 PetscSF sfMulti, sfInverse, sfCandidates; 1066 PetscInt *remoteOffsets; 1067 1068 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1069 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 1070 ierr = PetscSectionCreate(comm, &candidateRemoteSection);CHKERRQ(ierr); 1071 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);CHKERRQ(ierr); 1072 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);CHKERRQ(ierr); 1073 ierr = PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);CHKERRQ(ierr); 1074 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 1075 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1076 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1077 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 1078 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 1079 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1080 1081 ierr = PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");CHKERRQ(ierr); 1082 ierr = PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr); 1083 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr); 1084 } 1085 /* 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 */ 1086 { 1087 PetscHashIJKLRemote faceTable; 1088 PetscInt idx, idx2; 1089 1090 ierr = PetscHashIJKLRemoteCreate(&faceTable);CHKERRQ(ierr); 1091 /* There is a section point for every leaf attached to a given root point */ 1092 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1093 PetscInt deg; 1094 1095 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1096 PetscInt offset, dof, d; 1097 1098 ierr = PetscSectionGetDof(candidateRemoteSection, idx, &dof);CHKERRQ(ierr); 1099 ierr = PetscSectionGetOffset(candidateRemoteSection, idx, &offset);CHKERRQ(ierr); 1100 /* dof may include many faces from the remote process */ 1101 for (d = 0; d < dof; ++d) { 1102 const PetscInt hidx = offset+d; 1103 const PetscInt Np = candidatesRemote[hidx].index+1; 1104 const PetscSFNode rface = candidatesRemote[hidx+1]; 1105 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1106 PetscSFNode fcp0; 1107 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1108 const PetscInt *join = NULL; 1109 PetscHashIJKLRemoteKey key; 1110 PetscHashIter iter; 1111 PetscBool missing; 1112 PetscInt points[1024], p, joinSize; 1113 1114 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);} 1115 if (Np > 4) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%D, %D) at (%D, %D, %D) with %D cone points", rface.rank, rface.index, r, idx, d, Np); 1116 fcp0.rank = rank; 1117 fcp0.index = r; 1118 d += Np; 1119 /* Put remote face in hash table */ 1120 key.i = fcp0; 1121 key.j = fcone[0]; 1122 key.k = Np > 2 ? fcone[1] : pmax; 1123 key.l = Np > 3 ? fcone[2] : pmax; 1124 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1125 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1126 if (missing) { 1127 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1128 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1129 } else { 1130 PetscSFNode oface; 1131 1132 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1133 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1134 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1135 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1136 } 1137 } 1138 /* Check for local face */ 1139 points[0] = r; 1140 for (p = 1; p < Np; ++p) { 1141 ierr = DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]); 1142 if (ierr) break; /* We got a point not in our overlap */ 1143 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p]);CHKERRQ(ierr);} 1144 } 1145 if (ierr) continue; 1146 ierr = DMPlexGetJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1147 if (joinSize == 1) { 1148 PetscSFNode lface; 1149 PetscSFNode oface; 1150 1151 /* Always replace with local face */ 1152 lface.rank = rank; 1153 lface.index = join[0]; 1154 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1155 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);} 1156 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, lface);CHKERRQ(ierr); 1157 } 1158 ierr = DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1159 } 1160 } 1161 /* Put back faces for this root */ 1162 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1163 PetscInt offset, dof, d; 1164 1165 ierr = PetscSectionGetDof(candidateRemoteSection, idx2, &dof);CHKERRQ(ierr); 1166 ierr = PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);CHKERRQ(ierr); 1167 /* dof may include many faces from the remote process */ 1168 for (d = 0; d < dof; ++d) { 1169 const PetscInt hidx = offset+d; 1170 const PetscInt Np = candidatesRemote[hidx].index+1; 1171 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1172 PetscSFNode fcp0; 1173 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1174 PetscHashIJKLRemoteKey key; 1175 PetscHashIter iter; 1176 PetscBool missing; 1177 1178 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx);CHKERRQ(ierr);} 1179 if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np); 1180 fcp0.rank = rank; 1181 fcp0.index = r; 1182 d += Np; 1183 /* Find remote face in hash table */ 1184 key.i = fcp0; 1185 key.j = fcone[0]; 1186 key.k = Np > 2 ? fcone[1] : pmax; 1187 key.l = Np > 3 ? fcone[2] : pmax; 1188 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1189 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);} 1190 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1191 if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2); 1192 else {ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);CHKERRQ(ierr);} 1193 } 1194 } 1195 } 1196 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1197 ierr = PetscHashIJKLRemoteDestroy(&faceTable);CHKERRQ(ierr); 1198 } 1199 /* Step 4: Push back owned faces */ 1200 { 1201 PetscSF sfMulti, sfClaims, sfPointNew; 1202 PetscSFNode *remotePointsNew; 1203 PetscInt *remoteOffsets, *localPointsNew; 1204 PetscInt pStart, pEnd, r, NlNew, p; 1205 1206 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1207 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1208 ierr = PetscSectionCreate(comm, &claimSection);CHKERRQ(ierr); 1209 ierr = PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);CHKERRQ(ierr); 1210 ierr = PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1211 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr); 1212 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr); 1213 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1214 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1215 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1216 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1217 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1218 ierr = PetscObjectSetName((PetscObject) claimSection, "Claim Section");CHKERRQ(ierr); 1219 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr); 1220 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr); 1221 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1222 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1223 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1224 for (r = 0; r < Nr; ++r) { 1225 PetscInt dof, off, d; 1226 1227 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r);CHKERRQ(ierr);} 1228 ierr = PetscSectionGetDof(candidateSection, r, &dof);CHKERRQ(ierr); 1229 ierr = PetscSectionGetOffset(candidateSection, r, &off);CHKERRQ(ierr); 1230 for (d = 0; d < dof;) { 1231 if (claims[off+d].rank >= 0) { 1232 const PetscInt faceInd = off+d; 1233 const PetscInt Np = candidates[off+d].index; 1234 const PetscInt *join = NULL; 1235 PetscInt joinSize, points[1024], c; 1236 1237 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);CHKERRQ(ierr);} 1238 points[0] = r; 1239 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);} 1240 for (c = 0, d += 2; c < Np; ++c, ++d) { 1241 ierr = DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);CHKERRQ(ierr); 1242 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);} 1243 } 1244 ierr = DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1245 if (joinSize == 1) { 1246 if (claims[faceInd].rank == rank) { 1247 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0]);CHKERRQ(ierr);} 1248 } else { 1249 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);} 1250 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr); 1251 } 1252 } else { 1253 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank);CHKERRQ(ierr);} 1254 } 1255 ierr = DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1256 } else { 1257 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r);CHKERRQ(ierr);} 1258 d += claims[off+d].index+1; 1259 } 1260 } 1261 } 1262 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1263 /* Step 6) Create new pointSF from hashed claims */ 1264 ierr = PetscHMapIGetSize(claimshash, &NlNew);CHKERRQ(ierr); 1265 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1266 ierr = PetscMalloc1(Nl + NlNew, &localPointsNew);CHKERRQ(ierr); 1267 ierr = PetscMalloc1(Nl + NlNew, &remotePointsNew);CHKERRQ(ierr); 1268 for (l = 0; l < Nl; ++l) { 1269 localPointsNew[l] = localPoints[l]; 1270 remotePointsNew[l].index = remotePoints[l].index; 1271 remotePointsNew[l].rank = remotePoints[l].rank; 1272 } 1273 p = Nl; 1274 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1275 /* We sort new points, and assume they are numbered after all existing points */ 1276 ierr = PetscSortInt(NlNew, &localPointsNew[Nl]);CHKERRQ(ierr); 1277 for (p = Nl; p < Nl + NlNew; ++p) { 1278 PetscInt off; 1279 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &off);CHKERRQ(ierr); 1280 if (claims[off].rank < 0 || claims[off].index < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %D, (%D, %D)", localPointsNew[p], claims[off].rank, claims[off].index); 1281 remotePointsNew[p] = claims[off]; 1282 } 1283 ierr = PetscSFCreate(comm, &sfPointNew);CHKERRQ(ierr); 1284 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1285 ierr = PetscSFSetUp(sfPointNew);CHKERRQ(ierr); 1286 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1287 ierr = PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");CHKERRQ(ierr); 1288 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1289 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1290 } 1291 ierr = PetscHMapIJDestroy(&remoteHash);CHKERRQ(ierr); 1292 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1293 ierr = PetscSectionDestroy(&candidateRemoteSection);CHKERRQ(ierr); 1294 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1295 ierr = PetscFree(candidates);CHKERRQ(ierr); 1296 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1297 ierr = PetscFree(claims);CHKERRQ(ierr); 1298 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 /*@ 1303 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1304 1305 Collective on dm 1306 1307 Input Parameters: 1308 + dm - The DMPlex object with only cells and vertices 1309 - dmInt - The interpolated DM 1310 1311 Output Parameter: 1312 . dmInt - The complete DMPlex object 1313 1314 Level: intermediate 1315 1316 Notes: 1317 It does not copy over the coordinates. 1318 1319 Developer Notes: 1320 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL. 1321 1322 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1323 @*/ 1324 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1325 { 1326 DMPlexInterpolatedFlag interpolated; 1327 DM idm, odm = dm; 1328 PetscSF sfPoint; 1329 PetscInt depth, dim, d; 1330 const char *name; 1331 PetscBool flg=PETSC_TRUE; 1332 PetscErrorCode ierr; 1333 1334 PetscFunctionBegin; 1335 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1336 PetscValidPointer(dmInt, 2); 1337 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1338 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1339 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1340 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1341 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1342 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1343 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1344 idm = dm; 1345 } else { 1346 for (d = 1; d < dim; ++d) { 1347 /* Create interpolated mesh */ 1348 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1349 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1350 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1351 if (depth > 0) { 1352 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1353 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1354 { 1355 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1356 PetscInt nroots; 1357 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1358 if (nroots >= 0) {ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr);} 1359 } 1360 } 1361 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1362 odm = idm; 1363 } 1364 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1365 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1366 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1367 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr); 1368 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1369 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);} 1370 } 1371 { 1372 PetscBool isper; 1373 const PetscReal *maxCell, *L; 1374 const DMBoundaryType *bd; 1375 1376 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1377 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1378 } 1379 /* This function makes the mesh fully interpolated on all ranks */ 1380 { 1381 DM_Plex *plex = (DM_Plex *) idm->data; 1382 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1383 } 1384 *dmInt = idm; 1385 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1386 PetscFunctionReturn(0); 1387 } 1388 1389 /*@ 1390 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1391 1392 Collective on dmA 1393 1394 Input Parameter: 1395 . dmA - The DMPlex object with initial coordinates 1396 1397 Output Parameter: 1398 . dmB - The DMPlex object with copied coordinates 1399 1400 Level: intermediate 1401 1402 Note: This is typically used when adding pieces other than vertices to a mesh 1403 1404 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1405 @*/ 1406 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1407 { 1408 Vec coordinatesA, coordinatesB; 1409 VecType vtype; 1410 PetscSection coordSectionA, coordSectionB; 1411 PetscScalar *coordsA, *coordsB; 1412 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1413 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1414 PetscBool lc = PETSC_FALSE; 1415 PetscErrorCode ierr; 1416 1417 PetscFunctionBegin; 1418 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1419 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1420 if (dmA == dmB) PetscFunctionReturn(0); 1421 ierr = DMGetCoordinateDim(dmA, &cdim);CHKERRQ(ierr); 1422 ierr = DMSetCoordinateDim(dmB, cdim);CHKERRQ(ierr); 1423 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1424 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1425 if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB); 1426 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1427 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1428 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1429 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1430 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1431 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1432 if (!Nf) PetscFunctionReturn(0); 1433 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1434 if (!coordSectionB) { 1435 PetscInt dim; 1436 1437 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1438 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1439 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1440 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1441 } 1442 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1443 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1444 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1445 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1446 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1447 if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cells in first DM %D != %D in the second DM", cEndA-cStartA, cEndB-cStartB); 1448 cS = cS - cStartA + cStartB; 1449 cE = vEndB; 1450 lc = PETSC_TRUE; 1451 } else { 1452 cS = vStartB; 1453 cE = vEndB; 1454 } 1455 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1456 for (v = vStartB; v < vEndB; ++v) { 1457 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1458 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1459 } 1460 if (lc) { /* localized coordinates */ 1461 PetscInt c; 1462 1463 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1464 PetscInt dof; 1465 1466 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1467 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1468 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1469 } 1470 } 1471 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1472 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1473 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1474 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1475 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1476 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1477 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1478 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1479 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1480 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1481 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1482 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1483 for (v = 0; v < vEndB-vStartB; ++v) { 1484 PetscInt offA, offB; 1485 1486 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1487 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1488 for (d = 0; d < spaceDim; ++d) { 1489 coordsB[offB+d] = coordsA[offA+d]; 1490 } 1491 } 1492 if (lc) { /* localized coordinates */ 1493 PetscInt c; 1494 1495 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1496 PetscInt dof, offA, offB; 1497 1498 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1499 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1500 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1501 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1502 } 1503 } 1504 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1505 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1506 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1507 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /*@ 1512 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1513 1514 Collective on dm 1515 1516 Input Parameter: 1517 . dm - The complete DMPlex object 1518 1519 Output Parameter: 1520 . dmUnint - The DMPlex object with only cells and vertices 1521 1522 Level: intermediate 1523 1524 Notes: 1525 It does not copy over the coordinates. 1526 1527 Developer Notes: 1528 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1529 1530 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1531 @*/ 1532 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1533 { 1534 DMPlexInterpolatedFlag interpolated; 1535 DM udm; 1536 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1537 PetscErrorCode ierr; 1538 1539 PetscFunctionBegin; 1540 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1541 PetscValidPointer(dmUnint, 2); 1542 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1543 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1544 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1545 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1546 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1547 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1548 *dmUnint = dm; 1549 PetscFunctionReturn(0); 1550 } 1551 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1552 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1553 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1554 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1555 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1556 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1557 for (c = cStart; c < cEnd; ++c) { 1558 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1559 1560 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1561 for (cl = 0; cl < closureSize*2; cl += 2) { 1562 const PetscInt p = closure[cl]; 1563 1564 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1565 } 1566 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1567 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1568 maxConeSize = PetscMax(maxConeSize, coneSize); 1569 } 1570 ierr = DMSetUp(udm);CHKERRQ(ierr); 1571 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1572 for (c = cStart; c < cEnd; ++c) { 1573 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1574 1575 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1576 for (cl = 0; cl < closureSize*2; cl += 2) { 1577 const PetscInt p = closure[cl]; 1578 1579 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1580 } 1581 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1582 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1583 } 1584 ierr = PetscFree(cone);CHKERRQ(ierr); 1585 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1586 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1587 /* Reduce SF */ 1588 { 1589 PetscSF sfPoint, sfPointUn; 1590 const PetscSFNode *remotePoints; 1591 const PetscInt *localPoints; 1592 PetscSFNode *remotePointsUn; 1593 PetscInt *localPointsUn; 1594 PetscInt vEnd, numRoots, numLeaves, l; 1595 PetscInt numLeavesUn = 0, n = 0; 1596 PetscErrorCode ierr; 1597 1598 /* Get original SF information */ 1599 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1600 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1601 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1602 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1603 /* Allocate space for cells and vertices */ 1604 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1605 /* Fill in leaves */ 1606 if (vEnd >= 0) { 1607 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1608 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1609 for (l = 0; l < numLeaves; l++) { 1610 if (localPoints[l] < vEnd) { 1611 localPointsUn[n] = localPoints[l]; 1612 remotePointsUn[n].rank = remotePoints[l].rank; 1613 remotePointsUn[n].index = remotePoints[l].index; 1614 ++n; 1615 } 1616 } 1617 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1618 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1619 } 1620 } 1621 { 1622 PetscBool isper; 1623 const PetscReal *maxCell, *L; 1624 const DMBoundaryType *bd; 1625 1626 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1627 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1628 } 1629 /* This function makes the mesh fully uninterpolated on all ranks */ 1630 { 1631 DM_Plex *plex = (DM_Plex *) udm->data; 1632 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1633 } 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 if (flg != plex->interpolated) SETERRQ2(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);CHKERRQ(ierr); 1784 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRQ(ierr); 1785 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1786 if (debug) { 1787 PetscMPIInt rank; 1788 1789 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(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