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