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 /* 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 mainCone[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);CHKERRMPI(ierr); 644 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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,MPI_REPLACE);CHKERRQ(ierr); 674 ierr = PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks,MPI_REPLACE);CHKERRQ(ierr); 675 ierr = PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves,MPI_REPLACE);CHKERRQ(ierr); 676 ierr = PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks,MPI_REPLACE);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; mainCone 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 mainCone[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 mainCone[c] = rmine1[o+ind0];CHKERRQ(ierr); 704 } 705 if (debug) {ierr = PetscSynchronizedPrintf(comm, " mainCone=[%4D %4D]\n", mainCone[0], mainCone[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, mainCone);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);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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 static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared) 814 { 815 PetscSF sf; 816 const PetscInt *locals, *rootdegree; 817 PetscInt Nl, idx; 818 PetscErrorCode ierr; 819 820 PetscFunctionBegin; 821 *isShared = PETSC_FALSE; 822 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 823 ierr = PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL);CHKERRQ(ierr); 824 if (Nl < 0) PetscFunctionReturn(0); 825 ierr = PetscFindInt(p, Nl, locals, &idx);CHKERRQ(ierr); 826 if (idx >= 0) {*isShared = PETSC_TRUE; PetscFunctionReturn(0);} 827 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 828 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 829 if (rootdegree[p] > 0) *isShared = PETSC_TRUE; 830 PetscFunctionReturn(0); 831 } 832 833 static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared) 834 { 835 const PetscInt *cone; 836 PetscInt coneSize, c; 837 PetscBool cShared = PETSC_TRUE; 838 PetscErrorCode ierr; 839 840 PetscFunctionBegin; 841 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 842 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 843 for (c = 0; c < coneSize; ++c) { 844 PetscBool pointShared; 845 846 ierr = DMPlexPointIsShared(dm, cone[c], &pointShared);CHKERRQ(ierr); 847 cShared = (PetscBool) (cShared && pointShared); 848 } 849 *isShared = coneSize ? cShared : PETSC_FALSE; 850 PetscFunctionReturn(0); 851 } 852 853 static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin) 854 { 855 const PetscInt *cone; 856 PetscInt coneSize, c; 857 PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1}; 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 862 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 863 for (c = 0; c < coneSize; ++c) { 864 PetscSFNode rcp; 865 866 ierr = DMPlexMapToGlobalPoint(dm, cone[c], &rcp); 867 if (ierr) { 868 cmin = missing; 869 } else { 870 cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin; 871 } 872 } 873 *cpmin = coneSize ? cmin : missing; 874 PetscFunctionReturn(0); 875 } 876 877 /* 878 Each shared face has an entry in the candidates array: 879 (-1, coneSize-1), {(global cone point)} 880 where the set is missing the point p which we use as the key for the face 881 */ 882 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug) 883 { 884 MPI_Comm comm; 885 const PetscInt *support; 886 PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height; 887 PetscMPIInt rank; 888 PetscErrorCode ierr; 889 890 PetscFunctionBegin; 891 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 892 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 893 ierr = DMPlexGetOverlap(dm, &overlap);CHKERRQ(ierr); 894 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 895 ierr = DMPlexGetPointHeight(dm, p, &height);CHKERRQ(ierr); 896 if (!overlap && height <= cellHeight+1) { 897 /* cells can't be shared for non-overlapping meshes */ 898 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);} 899 PetscFunctionReturn(0); 900 } 901 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 902 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 903 if (candidates) {ierr = PetscSectionGetOffset(candidateSection, p, &off);CHKERRQ(ierr);} 904 for (s = 0; s < supportSize; ++s) { 905 const PetscInt face = support[s]; 906 const PetscInt *cone; 907 PetscSFNode cpmin={-1,-1}, rp={-1,-1}; 908 PetscInt coneSize, c, f; 909 PetscBool isShared = PETSC_FALSE; 910 PetscHashIJKey key; 911 912 /* Only add point once */ 913 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face);CHKERRQ(ierr);} 914 key.i = p; 915 key.j = face; 916 ierr = PetscHMapIJGet(faceHash, key, &f);CHKERRQ(ierr); 917 if (f >= 0) continue; 918 ierr = DMPlexConeIsShared(dm, face, &isShared);CHKERRQ(ierr); 919 ierr = DMPlexGetConeMinimum(dm, face, &cpmin);CHKERRQ(ierr); 920 ierr = DMPlexMapToGlobalPoint(dm, p, &rp);CHKERRQ(ierr); 921 if (debug) { 922 ierr = PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared);CHKERRQ(ierr); 923 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); 924 } 925 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 926 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 927 if (candidates) { 928 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank);CHKERRQ(ierr);} 929 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 930 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 931 candidates[off+idx].rank = -1; 932 candidates[off+idx++].index = coneSize-1; 933 candidates[off+idx].rank = rank; 934 candidates[off+idx++].index = face; 935 for (c = 0; c < coneSize; ++c) { 936 const PetscInt cp = cone[c]; 937 938 if (cp == p) continue; 939 ierr = DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);CHKERRQ(ierr); 940 if (debug) {ierr = PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);CHKERRQ(ierr);} 941 ++idx; 942 } 943 if (debug) {ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr);} 944 } else { 945 /* Add cone size to section */ 946 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face);CHKERRQ(ierr);} 947 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 948 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 949 ierr = PetscSectionAddDof(candidateSection, p, coneSize+1);CHKERRQ(ierr); 950 } 951 } 952 } 953 PetscFunctionReturn(0); 954 } 955 956 /*@ 957 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation 958 959 Collective on dm 960 961 Input Parameters: 962 + dm - The interpolated DM 963 - pointSF - The initial SF without interpolated points 964 965 Output Parameter: 966 . pointSF - The SF including interpolated points 967 968 Level: developer 969 970 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 971 972 .seealso: DMPlexInterpolate(), DMPlexUninterpolate() 973 @*/ 974 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 975 { 976 MPI_Comm comm; 977 PetscHMapIJ remoteHash; 978 PetscHMapI claimshash; 979 PetscSection candidateSection, candidateRemoteSection, claimSection; 980 PetscSFNode *candidates, *candidatesRemote, *claims; 981 const PetscInt *localPoints, *rootdegree; 982 const PetscSFNode *remotePoints; 983 PetscInt ov, Nr, r, Nl, l; 984 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 985 PetscBool flg, debug = PETSC_FALSE; 986 PetscMPIInt rank; 987 PetscErrorCode ierr; 988 989 PetscFunctionBegin; 990 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 991 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 2); 992 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 993 if (!flg) PetscFunctionReturn(0); 994 /* Set initial SF so that lower level queries work */ 995 ierr = DMSetPointSF(dm, pointSF);CHKERRQ(ierr); 996 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 997 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 998 ierr = DMPlexGetOverlap(dm, &ov);CHKERRQ(ierr); 999 if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1000 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr); 1001 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr); 1002 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr); 1003 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1004 /* Step 0: Precalculations */ 1005 ierr = PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);CHKERRQ(ierr); 1006 if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1007 ierr = PetscHMapIJCreate(&remoteHash);CHKERRQ(ierr); 1008 for (l = 0; l < Nl; ++l) { 1009 PetscHashIJKey key; 1010 key.i = remotePoints[l].index; 1011 key.j = remotePoints[l].rank; 1012 ierr = PetscHMapIJSet(remoteHash, key, l);CHKERRQ(ierr); 1013 } 1014 /* Compute root degree to identify shared points */ 1015 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 1016 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 1017 ierr = IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);CHKERRQ(ierr); 1018 /* 1019 1) Loop over each leaf point $p$ at depth $d$ in the SF 1020 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1021 \begin{itemize} 1022 \item all cone points of $f$ are shared 1023 \item $p$ is the cone point with smallest canonical number 1024 \end{itemize} 1025 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1026 \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 1027 \item Send the root face from the root back to all leaf process 1028 \item Leaf processes add the shared face to the SF 1029 */ 1030 /* Step 1: Construct section+SFNode array 1031 The section has entries for all shared faces for which we have a leaf point in the cone 1032 The array holds candidate shared faces, each face is refered to by the leaf point */ 1033 ierr = PetscSectionCreate(comm, &candidateSection);CHKERRQ(ierr); 1034 ierr = PetscSectionSetChart(candidateSection, 0, Nr);CHKERRQ(ierr); 1035 { 1036 PetscHMapIJ faceHash; 1037 1038 ierr = PetscHMapIJCreate(&faceHash);CHKERRQ(ierr); 1039 for (l = 0; l < Nl; ++l) { 1040 const PetscInt p = localPoints[l]; 1041 1042 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1043 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);CHKERRQ(ierr); 1044 } 1045 ierr = PetscHMapIJClear(faceHash);CHKERRQ(ierr); 1046 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 1047 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 1048 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 1049 for (l = 0; l < Nl; ++l) { 1050 const PetscInt p = localPoints[l]; 1051 1052 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1053 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);CHKERRQ(ierr); 1054 } 1055 ierr = PetscHMapIJDestroy(&faceHash);CHKERRQ(ierr); 1056 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1057 } 1058 ierr = PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");CHKERRQ(ierr); 1059 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr); 1060 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr); 1061 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1062 /* Note that this section is indexed by offsets into leaves, not by point number */ 1063 { 1064 PetscSF sfMulti, sfInverse, sfCandidates; 1065 PetscInt *remoteOffsets; 1066 1067 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1068 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 1069 ierr = PetscSectionCreate(comm, &candidateRemoteSection);CHKERRQ(ierr); 1070 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);CHKERRQ(ierr); 1071 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);CHKERRQ(ierr); 1072 ierr = PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);CHKERRQ(ierr); 1073 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 1074 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);CHKERRQ(ierr); 1075 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote,MPI_REPLACE);CHKERRQ(ierr); 1076 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 1077 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 1078 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1079 1080 ierr = PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");CHKERRQ(ierr); 1081 ierr = PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr); 1082 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr); 1083 } 1084 /* 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 */ 1085 { 1086 PetscHashIJKLRemote faceTable; 1087 PetscInt idx, idx2; 1088 1089 ierr = PetscHashIJKLRemoteCreate(&faceTable);CHKERRQ(ierr); 1090 /* There is a section point for every leaf attached to a given root point */ 1091 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1092 PetscInt deg; 1093 1094 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1095 PetscInt offset, dof, d; 1096 1097 ierr = PetscSectionGetDof(candidateRemoteSection, idx, &dof);CHKERRQ(ierr); 1098 ierr = PetscSectionGetOffset(candidateRemoteSection, idx, &offset);CHKERRQ(ierr); 1099 /* dof may include many faces from the remote process */ 1100 for (d = 0; d < dof; ++d) { 1101 const PetscInt hidx = offset+d; 1102 const PetscInt Np = candidatesRemote[hidx].index+1; 1103 const PetscSFNode rface = candidatesRemote[hidx+1]; 1104 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1105 PetscSFNode fcp0; 1106 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1107 const PetscInt *join = NULL; 1108 PetscHashIJKLRemoteKey key; 1109 PetscHashIter iter; 1110 PetscBool missing; 1111 PetscInt points[1024], p, joinSize; 1112 1113 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);} 1114 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); 1115 fcp0.rank = rank; 1116 fcp0.index = r; 1117 d += Np; 1118 /* Put remote face in hash table */ 1119 key.i = fcp0; 1120 key.j = fcone[0]; 1121 key.k = Np > 2 ? fcone[1] : pmax; 1122 key.l = Np > 3 ? fcone[2] : pmax; 1123 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1124 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1125 if (missing) { 1126 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1127 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1128 } else { 1129 PetscSFNode oface; 1130 1131 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1132 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1133 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1134 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1135 } 1136 } 1137 /* Check for local face */ 1138 points[0] = r; 1139 for (p = 1; p < Np; ++p) { 1140 ierr = DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]); 1141 if (ierr) break; /* We got a point not in our overlap */ 1142 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p]);CHKERRQ(ierr);} 1143 } 1144 if (ierr) continue; 1145 ierr = DMPlexGetJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1146 if (joinSize == 1) { 1147 PetscSFNode lface; 1148 PetscSFNode oface; 1149 1150 /* Always replace with local face */ 1151 lface.rank = rank; 1152 lface.index = join[0]; 1153 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1154 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);} 1155 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, lface);CHKERRQ(ierr); 1156 } 1157 ierr = DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1158 } 1159 } 1160 /* Put back faces for this root */ 1161 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1162 PetscInt offset, dof, d; 1163 1164 ierr = PetscSectionGetDof(candidateRemoteSection, idx2, &dof);CHKERRQ(ierr); 1165 ierr = PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);CHKERRQ(ierr); 1166 /* dof may include many faces from the remote process */ 1167 for (d = 0; d < dof; ++d) { 1168 const PetscInt hidx = offset+d; 1169 const PetscInt Np = candidatesRemote[hidx].index+1; 1170 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1171 PetscSFNode fcp0; 1172 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1173 PetscHashIJKLRemoteKey key; 1174 PetscHashIter iter; 1175 PetscBool missing; 1176 1177 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx);CHKERRQ(ierr);} 1178 if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np); 1179 fcp0.rank = rank; 1180 fcp0.index = r; 1181 d += Np; 1182 /* Find remote face in hash table */ 1183 key.i = fcp0; 1184 key.j = fcone[0]; 1185 key.k = Np > 2 ? fcone[1] : pmax; 1186 key.l = Np > 3 ? fcone[2] : pmax; 1187 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1188 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);} 1189 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1190 if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an associated face", r, idx2); 1191 else {ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);CHKERRQ(ierr);} 1192 } 1193 } 1194 } 1195 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1196 ierr = PetscHashIJKLRemoteDestroy(&faceTable);CHKERRQ(ierr); 1197 } 1198 /* Step 4: Push back owned faces */ 1199 { 1200 PetscSF sfMulti, sfClaims, sfPointNew; 1201 PetscSFNode *remotePointsNew; 1202 PetscInt *remoteOffsets, *localPointsNew; 1203 PetscInt pStart, pEnd, r, NlNew, p; 1204 1205 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1206 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1207 ierr = PetscSectionCreate(comm, &claimSection);CHKERRQ(ierr); 1208 ierr = PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);CHKERRQ(ierr); 1209 ierr = PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1210 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr); 1211 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr); 1212 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1213 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);CHKERRQ(ierr); 1214 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims,MPI_REPLACE);CHKERRQ(ierr); 1215 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1216 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1217 ierr = PetscObjectSetName((PetscObject) claimSection, "Claim Section");CHKERRQ(ierr); 1218 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr); 1219 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr); 1220 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1221 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1222 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1223 for (r = 0; r < Nr; ++r) { 1224 PetscInt dof, off, d; 1225 1226 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r);CHKERRQ(ierr);} 1227 ierr = PetscSectionGetDof(candidateSection, r, &dof);CHKERRQ(ierr); 1228 ierr = PetscSectionGetOffset(candidateSection, r, &off);CHKERRQ(ierr); 1229 for (d = 0; d < dof;) { 1230 if (claims[off+d].rank >= 0) { 1231 const PetscInt faceInd = off+d; 1232 const PetscInt Np = candidates[off+d].index; 1233 const PetscInt *join = NULL; 1234 PetscInt joinSize, points[1024], c; 1235 1236 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);CHKERRQ(ierr);} 1237 points[0] = r; 1238 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);} 1239 for (c = 0, d += 2; c < Np; ++c, ++d) { 1240 ierr = DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);CHKERRQ(ierr); 1241 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);} 1242 } 1243 ierr = DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1244 if (joinSize == 1) { 1245 if (claims[faceInd].rank == rank) { 1246 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0]);CHKERRQ(ierr);} 1247 } else { 1248 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);} 1249 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr); 1250 } 1251 } else { 1252 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank);CHKERRQ(ierr);} 1253 } 1254 ierr = DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1255 } else { 1256 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r);CHKERRQ(ierr);} 1257 d += claims[off+d].index+1; 1258 } 1259 } 1260 } 1261 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1262 /* Step 6) Create new pointSF from hashed claims */ 1263 ierr = PetscHMapIGetSize(claimshash, &NlNew);CHKERRQ(ierr); 1264 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1265 ierr = PetscMalloc1(Nl + NlNew, &localPointsNew);CHKERRQ(ierr); 1266 ierr = PetscMalloc1(Nl + NlNew, &remotePointsNew);CHKERRQ(ierr); 1267 for (l = 0; l < Nl; ++l) { 1268 localPointsNew[l] = localPoints[l]; 1269 remotePointsNew[l].index = remotePoints[l].index; 1270 remotePointsNew[l].rank = remotePoints[l].rank; 1271 } 1272 p = Nl; 1273 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1274 /* We sort new points, and assume they are numbered after all existing points */ 1275 ierr = PetscSortInt(NlNew, &localPointsNew[Nl]);CHKERRQ(ierr); 1276 for (p = Nl; p < Nl + NlNew; ++p) { 1277 PetscInt off; 1278 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &off);CHKERRQ(ierr); 1279 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); 1280 remotePointsNew[p] = claims[off]; 1281 } 1282 ierr = PetscSFCreate(comm, &sfPointNew);CHKERRQ(ierr); 1283 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1284 ierr = PetscSFSetUp(sfPointNew);CHKERRQ(ierr); 1285 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1286 ierr = PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");CHKERRQ(ierr); 1287 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1288 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1289 } 1290 ierr = PetscHMapIJDestroy(&remoteHash);CHKERRQ(ierr); 1291 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1292 ierr = PetscSectionDestroy(&candidateRemoteSection);CHKERRQ(ierr); 1293 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1294 ierr = PetscFree(candidates);CHKERRQ(ierr); 1295 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1296 ierr = PetscFree(claims);CHKERRQ(ierr); 1297 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /*@ 1302 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1303 1304 Collective on dm 1305 1306 Input Parameters: 1307 + dm - The DMPlex object with only cells and vertices 1308 - dmInt - The interpolated DM 1309 1310 Output Parameter: 1311 . dmInt - The complete DMPlex object 1312 1313 Level: intermediate 1314 1315 Notes: 1316 It does not copy over the coordinates. 1317 1318 Developer Notes: 1319 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL. 1320 1321 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1322 @*/ 1323 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1324 { 1325 DMPlexInterpolatedFlag interpolated; 1326 DM idm, odm = dm; 1327 PetscSF sfPoint; 1328 PetscInt depth, dim, d; 1329 const char *name; 1330 PetscBool flg=PETSC_TRUE; 1331 PetscErrorCode ierr; 1332 1333 PetscFunctionBegin; 1334 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1335 PetscValidPointer(dmInt, 2); 1336 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1337 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1338 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1339 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1340 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1341 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1342 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1343 idm = dm; 1344 } else { 1345 for (d = 1; d < dim; ++d) { 1346 /* Create interpolated mesh */ 1347 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1348 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1349 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1350 if (depth > 0) { 1351 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1352 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1353 { 1354 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1355 PetscInt nroots; 1356 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1357 if (nroots >= 0) {ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr);} 1358 } 1359 } 1360 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1361 odm = idm; 1362 } 1363 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1364 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1365 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1366 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr); 1367 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1368 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);} 1369 } 1370 { 1371 PetscBool isper; 1372 const PetscReal *maxCell, *L; 1373 const DMBoundaryType *bd; 1374 1375 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1376 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1377 } 1378 /* This function makes the mesh fully interpolated on all ranks */ 1379 { 1380 DM_Plex *plex = (DM_Plex *) idm->data; 1381 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1382 } 1383 *dmInt = idm; 1384 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1385 PetscFunctionReturn(0); 1386 } 1387 1388 /*@ 1389 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1390 1391 Collective on dmA 1392 1393 Input Parameter: 1394 . dmA - The DMPlex object with initial coordinates 1395 1396 Output Parameter: 1397 . dmB - The DMPlex object with copied coordinates 1398 1399 Level: intermediate 1400 1401 Note: This is typically used when adding pieces other than vertices to a mesh 1402 1403 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1404 @*/ 1405 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1406 { 1407 Vec coordinatesA, coordinatesB; 1408 VecType vtype; 1409 PetscSection coordSectionA, coordSectionB; 1410 PetscScalar *coordsA, *coordsB; 1411 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1412 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1413 PetscBool lc = PETSC_FALSE; 1414 PetscErrorCode ierr; 1415 1416 PetscFunctionBegin; 1417 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1418 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1419 if (dmA == dmB) PetscFunctionReturn(0); 1420 ierr = DMGetCoordinateDim(dmA, &cdim);CHKERRQ(ierr); 1421 ierr = DMSetCoordinateDim(dmB, cdim);CHKERRQ(ierr); 1422 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1423 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1424 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); 1425 /* Copy over discretization if it exists */ 1426 { 1427 DM cdmA, cdmB; 1428 PetscDS dsA, dsB; 1429 PetscObject objA, objB; 1430 PetscClassId idA, idB; 1431 const PetscScalar *constants; 1432 PetscInt cdim, Nc; 1433 1434 ierr = DMGetCoordinateDM(dmA, &cdmA);CHKERRQ(ierr); 1435 ierr = DMGetCoordinateDM(dmB, &cdmB);CHKERRQ(ierr); 1436 ierr = DMGetField(cdmA, 0, NULL, &objA);CHKERRQ(ierr); 1437 ierr = DMGetField(cdmB, 0, NULL, &objB);CHKERRQ(ierr); 1438 ierr = PetscObjectGetClassId(objA, &idA);CHKERRQ(ierr); 1439 ierr = PetscObjectGetClassId(objB, &idB);CHKERRQ(ierr); 1440 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1441 ierr = DMSetField(cdmB, 0, NULL, objA);CHKERRQ(ierr); 1442 ierr = DMCreateDS(cdmB);CHKERRQ(ierr); 1443 ierr = DMGetDS(cdmA, &dsA);CHKERRQ(ierr); 1444 ierr = DMGetDS(cdmB, &dsB);CHKERRQ(ierr); 1445 ierr = PetscDSGetCoordinateDimension(dsA, &cdim);CHKERRQ(ierr); 1446 ierr = PetscDSSetCoordinateDimension(dsB, cdim);CHKERRQ(ierr); 1447 ierr = PetscDSGetConstants(dsA, &Nc, &constants);CHKERRQ(ierr); 1448 ierr = PetscDSSetConstants(dsB, Nc, (PetscScalar *) constants);CHKERRQ(ierr); 1449 } 1450 } 1451 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1452 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1453 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1454 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1455 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1456 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1457 if (!Nf) PetscFunctionReturn(0); 1458 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1459 if (!coordSectionB) { 1460 PetscInt dim; 1461 1462 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1463 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1464 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1465 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1466 } 1467 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1468 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1469 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1470 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1471 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1472 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); 1473 cS = cS - cStartA + cStartB; 1474 cE = vEndB; 1475 lc = PETSC_TRUE; 1476 } else { 1477 cS = vStartB; 1478 cE = vEndB; 1479 } 1480 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1481 for (v = vStartB; v < vEndB; ++v) { 1482 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1483 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1484 } 1485 if (lc) { /* localized coordinates */ 1486 PetscInt c; 1487 1488 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1489 PetscInt dof; 1490 1491 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1492 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1493 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1494 } 1495 } 1496 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1497 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1498 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1499 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1500 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1501 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1502 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1503 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1504 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1505 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1506 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1507 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1508 for (v = 0; v < vEndB-vStartB; ++v) { 1509 PetscInt offA, offB; 1510 1511 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1512 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1513 for (d = 0; d < spaceDim; ++d) { 1514 coordsB[offB+d] = coordsA[offA+d]; 1515 } 1516 } 1517 if (lc) { /* localized coordinates */ 1518 PetscInt c; 1519 1520 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1521 PetscInt dof, offA, offB; 1522 1523 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1524 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1525 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1526 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1527 } 1528 } 1529 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1530 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1531 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1532 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1533 PetscFunctionReturn(0); 1534 } 1535 1536 /*@ 1537 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1538 1539 Collective on dm 1540 1541 Input Parameter: 1542 . dm - The complete DMPlex object 1543 1544 Output Parameter: 1545 . dmUnint - The DMPlex object with only cells and vertices 1546 1547 Level: intermediate 1548 1549 Notes: 1550 It does not copy over the coordinates. 1551 1552 Developer Notes: 1553 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1554 1555 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellListPetsc(), DMPlexCopyCoordinates() 1556 @*/ 1557 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1558 { 1559 DMPlexInterpolatedFlag interpolated; 1560 DM udm; 1561 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1562 PetscErrorCode ierr; 1563 1564 PetscFunctionBegin; 1565 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1566 PetscValidPointer(dmUnint, 2); 1567 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1568 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1569 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1570 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1571 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1572 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1573 *dmUnint = dm; 1574 PetscFunctionReturn(0); 1575 } 1576 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1577 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1578 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1579 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1580 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1581 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1582 for (c = cStart; c < cEnd; ++c) { 1583 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1584 1585 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1586 for (cl = 0; cl < closureSize*2; cl += 2) { 1587 const PetscInt p = closure[cl]; 1588 1589 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1590 } 1591 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1592 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1593 maxConeSize = PetscMax(maxConeSize, coneSize); 1594 } 1595 ierr = DMSetUp(udm);CHKERRQ(ierr); 1596 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1597 for (c = cStart; c < cEnd; ++c) { 1598 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1599 1600 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1601 for (cl = 0; cl < closureSize*2; cl += 2) { 1602 const PetscInt p = closure[cl]; 1603 1604 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1605 } 1606 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1607 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1608 } 1609 ierr = PetscFree(cone);CHKERRQ(ierr); 1610 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1611 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1612 /* Reduce SF */ 1613 { 1614 PetscSF sfPoint, sfPointUn; 1615 const PetscSFNode *remotePoints; 1616 const PetscInt *localPoints; 1617 PetscSFNode *remotePointsUn; 1618 PetscInt *localPointsUn; 1619 PetscInt vEnd, numRoots, numLeaves, l; 1620 PetscInt numLeavesUn = 0, n = 0; 1621 PetscErrorCode ierr; 1622 1623 /* Get original SF information */ 1624 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1625 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1626 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1627 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1628 /* Allocate space for cells and vertices */ 1629 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1630 /* Fill in leaves */ 1631 if (vEnd >= 0) { 1632 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1633 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1634 for (l = 0; l < numLeaves; l++) { 1635 if (localPoints[l] < vEnd) { 1636 localPointsUn[n] = localPoints[l]; 1637 remotePointsUn[n].rank = remotePoints[l].rank; 1638 remotePointsUn[n].index = remotePoints[l].index; 1639 ++n; 1640 } 1641 } 1642 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1643 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1644 } 1645 } 1646 { 1647 PetscBool isper; 1648 const PetscReal *maxCell, *L; 1649 const DMBoundaryType *bd; 1650 1651 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1652 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1653 } 1654 /* This function makes the mesh fully uninterpolated on all ranks */ 1655 { 1656 DM_Plex *plex = (DM_Plex *) udm->data; 1657 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1658 } 1659 *dmUnint = udm; 1660 PetscFunctionReturn(0); 1661 } 1662 1663 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1664 { 1665 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1666 MPI_Comm comm; 1667 PetscErrorCode ierr; 1668 1669 PetscFunctionBegin; 1670 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1671 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1672 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1673 1674 if (depth == dim) { 1675 *interpolated = DMPLEX_INTERPOLATED_FULL; 1676 if (!dim) goto finish; 1677 1678 /* Check points at height = dim are vertices (have no cones) */ 1679 ierr = DMPlexGetHeightStratum(dm, dim, &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 /* Check points at height < dim have cones */ 1689 for (h=0; h<dim; h++) { 1690 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 1691 for (p=pStart; p<pEnd; p++) { 1692 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1693 if (!coneSize) { 1694 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1695 goto finish; 1696 } 1697 } 1698 } 1699 } else if (depth == 1) { 1700 *interpolated = DMPLEX_INTERPOLATED_NONE; 1701 } else { 1702 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1703 } 1704 finish: 1705 PetscFunctionReturn(0); 1706 } 1707 1708 /*@ 1709 DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated. 1710 1711 Not Collective 1712 1713 Input Parameter: 1714 . dm - The DM object 1715 1716 Output Parameter: 1717 . interpolated - Flag whether the DM is interpolated 1718 1719 Level: intermediate 1720 1721 Notes: 1722 Unlike DMPlexIsInterpolatedCollective(), this is NOT collective 1723 so the results can be different on different ranks in special cases. 1724 However, DMPlexInterpolate() guarantees the result is the same on all. 1725 1726 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED. 1727 1728 Developer Notes: 1729 Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID. 1730 1731 If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called. 1732 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1733 DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL. 1734 1735 If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated. 1736 1737 DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL, 1738 and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1739 1740 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective() 1741 @*/ 1742 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1743 { 1744 DM_Plex *plex = (DM_Plex *) dm->data; 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1749 PetscValidPointer(interpolated,2); 1750 if (plex->interpolated < 0) { 1751 ierr = DMPlexIsInterpolated_Internal(dm, &plex->interpolated);CHKERRQ(ierr); 1752 } else if (PetscDefined (USE_DEBUG)) { 1753 DMPlexInterpolatedFlag flg; 1754 1755 ierr = DMPlexIsInterpolated_Internal(dm, &flg);CHKERRQ(ierr); 1756 if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1757 } 1758 *interpolated = plex->interpolated; 1759 PetscFunctionReturn(0); 1760 } 1761 1762 /*@ 1763 DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner). 1764 1765 Collective 1766 1767 Input Parameter: 1768 . dm - The DM object 1769 1770 Output Parameter: 1771 . interpolated - Flag whether the DM is interpolated 1772 1773 Level: intermediate 1774 1775 Notes: 1776 Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks. 1777 1778 This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks. 1779 1780 Developer Notes: 1781 Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID. 1782 1783 If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated. 1784 MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1785 if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED, 1786 otherwise sets plex->interpolatedCollective = plex->interpolated. 1787 1788 If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective. 1789 1790 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated() 1791 @*/ 1792 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1793 { 1794 DM_Plex *plex = (DM_Plex *) dm->data; 1795 PetscBool debug=PETSC_FALSE; 1796 PetscErrorCode ierr; 1797 1798 PetscFunctionBegin; 1799 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1800 PetscValidPointer(interpolated,2); 1801 ierr = PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);CHKERRQ(ierr); 1802 if (plex->interpolatedCollective < 0) { 1803 DMPlexInterpolatedFlag min, max; 1804 MPI_Comm comm; 1805 1806 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1807 ierr = DMPlexIsInterpolated(dm, &plex->interpolatedCollective);CHKERRQ(ierr); 1808 ierr = MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);CHKERRMPI(ierr); 1809 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRMPI(ierr); 1810 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1811 if (debug) { 1812 PetscMPIInt rank; 1813 1814 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1815 ierr = PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);CHKERRQ(ierr); 1816 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 1817 } 1818 } 1819 *interpolated = plex->interpolatedCollective; 1820 PetscFunctionReturn(0); 1821 } 1822