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