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