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 static PetscErrorCode DMPlexCheckSharedFaces(DM dm) 1039 { 1040 PetscInt fStart, fEnd, f; 1041 PetscErrorCode ierr; 1042 1043 PetscFunctionBegin; 1044 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1045 for (f = fStart; f < fEnd; ++f) { 1046 const PetscInt *cone; 1047 PetscInt coneSize, c; 1048 PetscBool faceShared, shouldShare = PETSC_TRUE; 1049 1050 ierr = DMPlexGetConeSize(dm, f, &coneSize);CHKERRQ(ierr); 1051 ierr = DMPlexGetCone(dm, f, &cone);CHKERRQ(ierr); 1052 for (c = 0; c < coneSize; ++c) { 1053 PetscBool pointShared; 1054 1055 ierr = DMPlexPointIsShared(dm, cone[c], &pointShared);CHKERRQ(ierr); 1056 shouldShare = (PetscBool) (shouldShare && pointShared); 1057 } 1058 ierr = DMPlexPointIsShared(dm, f, &faceShared);CHKERRQ(ierr); 1059 if (faceShared && !shouldShare) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Face %D is shared, but its cone is not completely shared", f); 1060 if (!faceShared && shouldShare) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Face %D is not shared, but its cone is completely shared", f); 1061 } 1062 PetscFunctionReturn(0); 1063 } 1064 1065 /* 1066 Each shared face has an entry in the candidates array: 1067 (-1, coneSize-1), {(global cone point)} 1068 where the set is missing the point p which we use as the key for the face 1069 */ 1070 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug) 1071 { 1072 MPI_Comm comm; 1073 const PetscInt *support; 1074 PetscInt supportSize, s, off = 0, idx = 0; 1075 PetscMPIInt rank; 1076 PetscErrorCode ierr; 1077 1078 PetscFunctionBegin; 1079 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 1080 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1081 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 1082 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 1083 if (candidates) {ierr = PetscSectionGetOffset(candidateSection, p, &off);CHKERRQ(ierr);} 1084 for (s = 0; s < supportSize; ++s) { 1085 const PetscInt face = support[s]; 1086 const PetscInt *cone; 1087 PetscSFNode cpmin={-1,-1}, rp={-1,-1}; 1088 PetscInt coneSize, c, f; 1089 PetscBool isShared = PETSC_FALSE; 1090 PetscHashIJKey key; 1091 1092 /* Only add point once */ 1093 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Support face %D\n", rank, face);CHKERRQ(ierr);} 1094 key.i = p; 1095 key.j = face; 1096 ierr = PetscHMapIJGet(faceHash, key, &f);CHKERRQ(ierr); 1097 if (f >= 0) continue; 1098 ierr = DMPlexConeIsShared(dm, face, &isShared);CHKERRQ(ierr); 1099 ierr = DMPlexGetConeMinimum(dm, face, &cpmin);CHKERRQ(ierr); 1100 ierr = DMPlexMapToGlobalPoint(dm, p, &rp);CHKERRQ(ierr); 1101 if (debug) { 1102 ierr = PetscSynchronizedPrintf(comm, "[%d] Face point %D is shared: %d\n", rank, face, (int) isShared);CHKERRQ(ierr); 1103 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); 1104 } 1105 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 1106 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 1107 if (candidates) { 1108 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Adding shared face %D at idx %D\n[%d] ", rank, face, idx, rank);CHKERRQ(ierr);} 1109 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1110 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 1111 candidates[off+idx].rank = -1; 1112 candidates[off+idx++].index = coneSize-1; 1113 candidates[off+idx].rank = rank; 1114 candidates[off+idx++].index = face; 1115 for (c = 0; c < coneSize; ++c) { 1116 const PetscInt cp = cone[c]; 1117 1118 if (cp == p) continue; 1119 ierr = DMPlexMapToGlobalPoint(dm, cp, &candidates[off+idx]);CHKERRQ(ierr); 1120 if (debug) {ierr = PetscSynchronizedPrintf(comm, " (%D,%D)", candidates[off+idx].rank, candidates[off+idx].index);CHKERRQ(ierr);} 1121 ++idx; 1122 } 1123 if (debug) {ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr);} 1124 } else { 1125 /* Add cone size to section */ 1126 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %D\n", rank, face);CHKERRQ(ierr);} 1127 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1128 ierr = PetscHMapIJSet(faceHash, key, p);CHKERRQ(ierr); 1129 ierr = PetscSectionAddDof(candidateSection, p, coneSize+1);CHKERRQ(ierr); 1130 } 1131 } 1132 } 1133 PetscFunctionReturn(0); 1134 } 1135 1136 /*@ 1137 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation 1138 1139 Collective on dm 1140 1141 Input Parameters: 1142 + dm - The interpolated DM 1143 - pointSF - The initial SF without interpolated points 1144 1145 Output Parameter: 1146 . pointSF - The SF including interpolated points 1147 1148 Level: developer 1149 1150 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 1151 1152 .seealso: DMPlexInterpolate(), DMPlexUninterpolate() 1153 @*/ 1154 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 1155 { 1156 MPI_Comm comm; 1157 PetscHMapIJ remoteHash; 1158 PetscHMapI claimshash; 1159 PetscSection candidateSection, candidateRemoteSection, claimSection; 1160 PetscSFNode *candidates, *candidatesRemote, *claims; 1161 const PetscInt *localPoints, *rootdegree; 1162 const PetscSFNode *remotePoints; 1163 PetscInt ov, Nr, r, Nl, l; 1164 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 1165 PetscBool flg, debug = PETSC_FALSE; 1166 PetscMPIInt rank; 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1171 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 3); 1172 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 1173 if (!flg) PetscFunctionReturn(0); 1174 /* Set initial SF so that lower level queries work */ 1175 ierr = DMSetPointSF(dm, pointSF);CHKERRQ(ierr); 1176 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 1177 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1178 ierr = DMPlexGetOverlap(dm, &ov);CHKERRQ(ierr); 1179 if (ov) SETERRQ(comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1180 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr); 1181 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr); 1182 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr); 1183 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1184 /* Step 0: Precalculations */ 1185 ierr = PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints);CHKERRQ(ierr); 1186 if (Nr < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1187 ierr = PetscHMapIJCreate(&remoteHash);CHKERRQ(ierr); 1188 for (l = 0; l < Nl; ++l) { 1189 PetscHashIJKey key; 1190 key.i = remotePoints[l].index; 1191 key.j = remotePoints[l].rank; 1192 ierr = PetscHMapIJSet(remoteHash, key, l);CHKERRQ(ierr); 1193 } 1194 /* Compute root degree to identify shared points */ 1195 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 1196 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 1197 ierr = IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree);CHKERRQ(ierr); 1198 /* 1199 1) Loop over each leaf point $p$ at depth $d$ in the SF 1200 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1201 \begin{itemize} 1202 \item all cone points of $f$ are shared 1203 \item $p$ is the cone point with smallest canonical number 1204 \end{itemize} 1205 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1206 \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 1207 \item Send the root face from the root back to all leaf process 1208 \item Leaf processes add the shared face to the SF 1209 */ 1210 /* Step 1: Construct section+SFNode array 1211 The section has entries for all shared faces for which we have a leaf point in the cone 1212 The array holds candidate shared faces, each face is refered to by the leaf point */ 1213 ierr = PetscSectionCreate(comm, &candidateSection);CHKERRQ(ierr); 1214 ierr = PetscSectionSetChart(candidateSection, 0, Nr);CHKERRQ(ierr); 1215 { 1216 PetscHMapIJ faceHash; 1217 1218 ierr = PetscHMapIJCreate(&faceHash);CHKERRQ(ierr); 1219 for (l = 0; l < Nl; ++l) { 1220 const PetscInt p = localPoints[l]; 1221 1222 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1223 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug);CHKERRQ(ierr); 1224 } 1225 ierr = PetscHMapIJClear(faceHash);CHKERRQ(ierr); 1226 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 1227 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 1228 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 1229 for (l = 0; l < Nl; ++l) { 1230 const PetscInt p = localPoints[l]; 1231 1232 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %D\n", rank, p);CHKERRQ(ierr);} 1233 ierr = DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug);CHKERRQ(ierr); 1234 } 1235 ierr = PetscHMapIJDestroy(&faceHash);CHKERRQ(ierr); 1236 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1237 } 1238 ierr = PetscObjectSetName((PetscObject) candidateSection, "Candidate Section");CHKERRQ(ierr); 1239 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr); 1240 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr); 1241 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1242 /* Note that this section is indexed by offsets into leaves, not by point number */ 1243 { 1244 PetscSF sfMulti, sfInverse, sfCandidates; 1245 PetscInt *remoteOffsets; 1246 1247 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1248 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 1249 ierr = PetscSectionCreate(comm, &candidateRemoteSection);CHKERRQ(ierr); 1250 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection);CHKERRQ(ierr); 1251 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates);CHKERRQ(ierr); 1252 ierr = PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize);CHKERRQ(ierr); 1253 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 1254 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1255 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1256 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 1257 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 1258 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1259 1260 ierr = PetscObjectSetName((PetscObject) candidateRemoteSection, "Remote Candidate Section");CHKERRQ(ierr); 1261 ierr = PetscObjectViewFromOptions((PetscObject) candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr); 1262 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr); 1263 } 1264 /* 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 */ 1265 { 1266 PetscHashIJKLRemote faceTable; 1267 PetscInt idx, idx2; 1268 1269 ierr = PetscHashIJKLRemoteCreate(&faceTable);CHKERRQ(ierr); 1270 /* There is a section point for every leaf attached to a given root point */ 1271 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1272 PetscInt deg; 1273 1274 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1275 PetscInt offset, dof, d; 1276 1277 ierr = PetscSectionGetDof(candidateRemoteSection, idx, &dof);CHKERRQ(ierr); 1278 ierr = PetscSectionGetOffset(candidateRemoteSection, idx, &offset);CHKERRQ(ierr); 1279 /* dof may include many faces from the remote process */ 1280 for (d = 0; d < dof; ++d) { 1281 const PetscInt hidx = offset+d; 1282 const PetscInt Np = candidatesRemote[hidx].index+1; 1283 const PetscSFNode rface = candidatesRemote[hidx+1]; 1284 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1285 PetscSFNode fcp0; 1286 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1287 const PetscInt *join = NULL; 1288 PetscHashIJKLRemoteKey key; 1289 PetscHashIter iter; 1290 PetscBool missing; 1291 PetscInt points[1024], p, joinSize; 1292 1293 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);} 1294 if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np); 1295 fcp0.rank = rank; 1296 fcp0.index = r; 1297 d += Np; 1298 /* Put remote face in hash table */ 1299 key.i = fcp0; 1300 key.j = fcone[0]; 1301 key.k = Np > 2 ? fcone[1] : pmax; 1302 key.l = Np > 3 ? fcone[2] : pmax; 1303 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1304 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1305 if (missing) { 1306 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Setting remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1307 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1308 } else { 1309 PetscSFNode oface; 1310 1311 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1312 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1313 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Replacing with remote face (%D, %D)\n", rank, rface.index, rface.rank);CHKERRQ(ierr);} 1314 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, rface);CHKERRQ(ierr); 1315 } 1316 } 1317 /* Check for local face */ 1318 points[0] = r; 1319 for (p = 1; p < Np; ++p) { 1320 ierr = DMPlexMapToLocalPoint(dm, remoteHash, fcone[p-1], &points[p]); 1321 if (ierr) break; /* We got a point not in our overlap */ 1322 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p]);CHKERRQ(ierr);} 1323 } 1324 if (ierr) continue; 1325 ierr = DMPlexGetJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1326 if (joinSize == 1) { 1327 PetscSFNode lface; 1328 PetscSFNode oface; 1329 1330 /* Always replace with local face */ 1331 lface.rank = rank; 1332 lface.index = join[0]; 1333 ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &oface);CHKERRQ(ierr); 1334 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);} 1335 ierr = PetscHashIJKLRemoteIterSet(faceTable, iter, lface);CHKERRQ(ierr); 1336 } 1337 ierr = DMPlexRestoreJoin(dm, Np, points, &joinSize, &join);CHKERRQ(ierr); 1338 } 1339 } 1340 /* Put back faces for this root */ 1341 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1342 PetscInt offset, dof, d; 1343 1344 ierr = PetscSectionGetDof(candidateRemoteSection, idx2, &dof);CHKERRQ(ierr); 1345 ierr = PetscSectionGetOffset(candidateRemoteSection, idx2, &offset);CHKERRQ(ierr); 1346 /* dof may include many faces from the remote process */ 1347 for (d = 0; d < dof; ++d) { 1348 const PetscInt hidx = offset+d; 1349 const PetscInt Np = candidatesRemote[hidx].index+1; 1350 const PetscSFNode *fcone = &candidatesRemote[hidx+2]; 1351 PetscSFNode fcp0; 1352 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1353 PetscHashIJKLRemoteKey key; 1354 PetscHashIter iter; 1355 PetscBool missing; 1356 1357 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Entering face at (%D, %D)\n", rank, r, idx);CHKERRQ(ierr);} 1358 if (Np > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %D cone points", Np); 1359 fcp0.rank = rank; 1360 fcp0.index = r; 1361 d += Np; 1362 /* Find remote face in hash table */ 1363 key.i = fcp0; 1364 key.j = fcone[0]; 1365 key.k = Np > 2 ? fcone[1] : pmax; 1366 key.l = Np > 3 ? fcone[2] : pmax; 1367 ierr = PetscSortSFNode(Np, (PetscSFNode *) &key);CHKERRQ(ierr); 1368 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);} 1369 ierr = PetscHashIJKLRemotePut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 1370 if (missing) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %D Idx %D ought to have an assoicated face", r, idx2); 1371 else {ierr = PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx]);CHKERRQ(ierr);} 1372 } 1373 } 1374 } 1375 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1376 ierr = PetscHashIJKLRemoteDestroy(&faceTable);CHKERRQ(ierr); 1377 } 1378 /* Step 4: Push back owned faces */ 1379 { 1380 PetscSF sfMulti, sfClaims, sfPointNew; 1381 PetscSFNode *remotePointsNew; 1382 PetscInt *remoteOffsets, *localPointsNew; 1383 PetscInt pStart, pEnd, r, NlNew, p; 1384 1385 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1386 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1387 ierr = PetscSectionCreate(comm, &claimSection);CHKERRQ(ierr); 1388 ierr = PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection);CHKERRQ(ierr); 1389 ierr = PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1390 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr); 1391 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr); 1392 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1393 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1394 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1395 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1396 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1397 ierr = PetscObjectSetName((PetscObject) claimSection, "Claim Section");CHKERRQ(ierr); 1398 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr); 1399 ierr = SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr); 1400 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1401 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1402 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1403 for (r = 0; r < Nr; ++r) { 1404 PetscInt dof, off, d; 1405 1406 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %D\n", rank, r);CHKERRQ(ierr);} 1407 ierr = PetscSectionGetDof(candidateSection, r, &dof);CHKERRQ(ierr); 1408 ierr = PetscSectionGetOffset(candidateSection, r, &off);CHKERRQ(ierr); 1409 for (d = 0; d < dof;) { 1410 if (claims[off+d].rank >= 0) { 1411 const PetscInt faceInd = off+d; 1412 const PetscInt Np = candidates[off+d].index; 1413 const PetscInt *join = NULL; 1414 PetscInt joinSize, points[1024], c; 1415 1416 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%D, %D)\n", rank, claims[faceInd].rank, claims[faceInd].index);CHKERRQ(ierr);} 1417 points[0] = r; 1418 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);} 1419 for (c = 0, d += 2; c < Np; ++c, ++d) { 1420 ierr = DMPlexMapToLocalPoint(dm, remoteHash, candidates[off+d], &points[c+1]);CHKERRQ(ierr); 1421 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);} 1422 } 1423 ierr = DMPlexGetJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1424 if (joinSize == 1) { 1425 if (claims[faceInd].rank == rank) { 1426 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %D for non-remote partner\n", rank, join[0]);CHKERRQ(ierr);} 1427 } else { 1428 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);} 1429 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr); 1430 } 1431 } else { 1432 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank);CHKERRQ(ierr);} 1433 } 1434 ierr = DMPlexRestoreJoin(dm, Np+1, points, &joinSize, &join);CHKERRQ(ierr); 1435 } else { 1436 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] No claim for point %D\n", rank, r);CHKERRQ(ierr);} 1437 d += claims[off+d].index+1; 1438 } 1439 } 1440 } 1441 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 1442 /* Step 6) Create new pointSF from hashed claims */ 1443 ierr = PetscHMapIGetSize(claimshash, &NlNew);CHKERRQ(ierr); 1444 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1445 ierr = PetscMalloc1(Nl + NlNew, &localPointsNew);CHKERRQ(ierr); 1446 ierr = PetscMalloc1(Nl + NlNew, &remotePointsNew);CHKERRQ(ierr); 1447 for (l = 0; l < Nl; ++l) { 1448 localPointsNew[l] = localPoints[l]; 1449 remotePointsNew[l].index = remotePoints[l].index; 1450 remotePointsNew[l].rank = remotePoints[l].rank; 1451 } 1452 p = Nl; 1453 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1454 /* We sort new points, and assume they are numbered after all existing points */ 1455 ierr = PetscSortInt(NlNew, &localPointsNew[Nl]);CHKERRQ(ierr); 1456 for (p = Nl; p < Nl + NlNew; ++p) { 1457 PetscInt off; 1458 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &off);CHKERRQ(ierr); 1459 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); 1460 remotePointsNew[p] = claims[off]; 1461 } 1462 ierr = PetscSFCreate(comm, &sfPointNew);CHKERRQ(ierr); 1463 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, Nl+NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1464 ierr = PetscSFSetUp(sfPointNew);CHKERRQ(ierr); 1465 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1466 ierr = PetscObjectViewFromOptions((PetscObject) sfPointNew, NULL, "-petscsf_interp_view");CHKERRQ(ierr); 1467 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1468 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1469 } 1470 ierr = PetscHMapIJDestroy(&remoteHash);CHKERRQ(ierr); 1471 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1472 ierr = PetscSectionDestroy(&candidateRemoteSection);CHKERRQ(ierr); 1473 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1474 ierr = PetscFree(candidates);CHKERRQ(ierr); 1475 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1476 ierr = PetscFree(claims);CHKERRQ(ierr); 1477 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 /*@ 1482 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1483 1484 Collective on dm 1485 1486 Input Parameters: 1487 + dm - The DMPlex object with only cells and vertices 1488 - dmInt - The interpolated DM 1489 1490 Output Parameter: 1491 . dmInt - The complete DMPlex object 1492 1493 Level: intermediate 1494 1495 Notes: 1496 It does not copy over the coordinates. 1497 1498 Developer Notes: 1499 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL. 1500 1501 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1502 @*/ 1503 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1504 { 1505 DMPlexInterpolatedFlag interpolated; 1506 DM idm, odm = dm; 1507 PetscSF sfPoint; 1508 PetscInt depth, dim, d; 1509 const char *name; 1510 PetscBool flg=PETSC_TRUE; 1511 PetscErrorCode ierr; 1512 1513 PetscFunctionBegin; 1514 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1515 PetscValidPointer(dmInt, 2); 1516 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1517 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1518 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1519 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1520 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1521 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1522 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1523 idm = dm; 1524 } else { 1525 for (d = 1; d < dim; ++d) { 1526 /* Create interpolated mesh */ 1527 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1528 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1529 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1530 if (depth > 0) { 1531 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1532 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1533 { 1534 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1535 PetscInt nroots; 1536 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1537 if (nroots >= 0) {ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr);} 1538 } 1539 } 1540 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1541 odm = idm; 1542 } 1543 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1544 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1545 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1546 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr); 1547 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1548 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);} 1549 } 1550 { 1551 PetscBool isper; 1552 const PetscReal *maxCell, *L; 1553 const DMBoundaryType *bd; 1554 1555 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1556 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1557 } 1558 /* This function makes the mesh fully interpolated on all ranks */ 1559 { 1560 DM_Plex *plex = (DM_Plex *) idm->data; 1561 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1562 } 1563 *dmInt = idm; 1564 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 /*@ 1569 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1570 1571 Collective on dmA 1572 1573 Input Parameter: 1574 . dmA - The DMPlex object with initial coordinates 1575 1576 Output Parameter: 1577 . dmB - The DMPlex object with copied coordinates 1578 1579 Level: intermediate 1580 1581 Note: This is typically used when adding pieces other than vertices to a mesh 1582 1583 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1584 @*/ 1585 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1586 { 1587 Vec coordinatesA, coordinatesB; 1588 VecType vtype; 1589 PetscSection coordSectionA, coordSectionB; 1590 PetscScalar *coordsA, *coordsB; 1591 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1592 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 1593 PetscBool lc = PETSC_FALSE; 1594 PetscErrorCode ierr; 1595 1596 PetscFunctionBegin; 1597 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1598 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1599 if (dmA == dmB) PetscFunctionReturn(0); 1600 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1601 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1602 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); 1603 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1604 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1605 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1606 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1607 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1608 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1609 if (!Nf) PetscFunctionReturn(0); 1610 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1611 if (!coordSectionB) { 1612 PetscInt dim; 1613 1614 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1615 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1616 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1617 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1618 } 1619 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1620 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1621 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1622 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1623 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1624 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); 1625 cS = cS - cStartA + cStartB; 1626 cE = vEndB; 1627 lc = PETSC_TRUE; 1628 } else { 1629 cS = vStartB; 1630 cE = vEndB; 1631 } 1632 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1633 for (v = vStartB; v < vEndB; ++v) { 1634 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1635 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1636 } 1637 if (lc) { /* localized coordinates */ 1638 PetscInt c; 1639 1640 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1641 PetscInt dof; 1642 1643 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1644 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1645 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1646 } 1647 } 1648 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1649 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1650 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1651 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1652 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1653 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1654 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1655 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1656 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1657 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1658 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1659 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1660 for (v = 0; v < vEndB-vStartB; ++v) { 1661 PetscInt offA, offB; 1662 1663 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1664 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1665 for (d = 0; d < spaceDim; ++d) { 1666 coordsB[offB+d] = coordsA[offA+d]; 1667 } 1668 } 1669 if (lc) { /* localized coordinates */ 1670 PetscInt c; 1671 1672 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1673 PetscInt dof, offA, offB; 1674 1675 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1676 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1677 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1678 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1679 } 1680 } 1681 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1682 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1683 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1684 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 /*@ 1689 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1690 1691 Collective on dm 1692 1693 Input Parameter: 1694 . dm - The complete DMPlex object 1695 1696 Output Parameter: 1697 . dmUnint - The DMPlex object with only cells and vertices 1698 1699 Level: intermediate 1700 1701 Notes: 1702 It does not copy over the coordinates. 1703 1704 Developer Notes: 1705 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1706 1707 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1708 @*/ 1709 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1710 { 1711 DMPlexInterpolatedFlag interpolated; 1712 DM udm; 1713 PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone; 1714 PetscErrorCode ierr; 1715 1716 PetscFunctionBegin; 1717 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1718 PetscValidPointer(dmUnint, 2); 1719 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1720 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1721 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1722 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1723 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1724 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1725 *dmUnint = dm; 1726 PetscFunctionReturn(0); 1727 } 1728 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1729 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1730 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1731 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1732 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1733 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1734 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1735 for (c = cStart; c < cEnd; ++c) { 1736 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1737 1738 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1739 for (cl = 0; cl < closureSize*2; cl += 2) { 1740 const PetscInt p = closure[cl]; 1741 1742 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1743 } 1744 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1745 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1746 maxConeSize = PetscMax(maxConeSize, coneSize); 1747 } 1748 ierr = DMSetUp(udm);CHKERRQ(ierr); 1749 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1750 for (c = cStart; c < cEnd; ++c) { 1751 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1752 1753 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1754 for (cl = 0; cl < closureSize*2; cl += 2) { 1755 const PetscInt p = closure[cl]; 1756 1757 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1758 } 1759 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1760 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1761 } 1762 ierr = PetscFree(cone);CHKERRQ(ierr); 1763 ierr = DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1764 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1765 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1766 /* Reduce SF */ 1767 { 1768 PetscSF sfPoint, sfPointUn; 1769 const PetscSFNode *remotePoints; 1770 const PetscInt *localPoints; 1771 PetscSFNode *remotePointsUn; 1772 PetscInt *localPointsUn; 1773 PetscInt vEnd, numRoots, numLeaves, l; 1774 PetscInt numLeavesUn = 0, n = 0; 1775 PetscErrorCode ierr; 1776 1777 /* Get original SF information */ 1778 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1779 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1780 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1781 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1782 /* Allocate space for cells and vertices */ 1783 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1784 /* Fill in leaves */ 1785 if (vEnd >= 0) { 1786 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1787 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1788 for (l = 0; l < numLeaves; l++) { 1789 if (localPoints[l] < vEnd) { 1790 localPointsUn[n] = localPoints[l]; 1791 remotePointsUn[n].rank = remotePoints[l].rank; 1792 remotePointsUn[n].index = remotePoints[l].index; 1793 ++n; 1794 } 1795 } 1796 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1797 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1798 } 1799 } 1800 { 1801 PetscBool isper; 1802 const PetscReal *maxCell, *L; 1803 const DMBoundaryType *bd; 1804 1805 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1806 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1807 } 1808 /* This function makes the mesh fully uninterpolated on all ranks */ 1809 { 1810 DM_Plex *plex = (DM_Plex *) udm->data; 1811 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1812 } 1813 *dmUnint = udm; 1814 PetscFunctionReturn(0); 1815 } 1816 1817 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1818 { 1819 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1820 MPI_Comm comm; 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1825 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1826 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1827 1828 if (depth == dim) { 1829 *interpolated = DMPLEX_INTERPOLATED_FULL; 1830 if (!dim) goto finish; 1831 1832 /* Check points at height = dim are vertices (have no cones) */ 1833 ierr = DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);CHKERRQ(ierr); 1834 for (p=pStart; p<pEnd; p++) { 1835 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1836 if (coneSize) { 1837 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1838 goto finish; 1839 } 1840 } 1841 1842 /* Check points at height < dim have cones */ 1843 for (h=0; h<dim; h++) { 1844 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 1845 for (p=pStart; p<pEnd; p++) { 1846 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1847 if (!coneSize) { 1848 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1849 goto finish; 1850 } 1851 } 1852 } 1853 } else if (depth == 1) { 1854 *interpolated = DMPLEX_INTERPOLATED_NONE; 1855 } else { 1856 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1857 } 1858 finish: 1859 PetscFunctionReturn(0); 1860 } 1861 1862 /*@ 1863 DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated. 1864 1865 Not Collective 1866 1867 Input Parameter: 1868 . dm - The DM object 1869 1870 Output Parameter: 1871 . interpolated - Flag whether the DM is interpolated 1872 1873 Level: intermediate 1874 1875 Notes: 1876 Unlike DMPlexIsInterpolatedCollective(), this is NOT collective 1877 so the results can be different on different ranks in special cases. 1878 However, DMPlexInterpolate() guarantees the result is the same on all. 1879 1880 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED. 1881 1882 Developer Notes: 1883 Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID. 1884 1885 If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called. 1886 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1887 DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL. 1888 1889 If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated. 1890 1891 DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL, 1892 and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1893 1894 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective() 1895 @*/ 1896 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1897 { 1898 DM_Plex *plex = (DM_Plex *) dm->data; 1899 PetscErrorCode ierr; 1900 1901 PetscFunctionBegin; 1902 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1903 PetscValidPointer(interpolated,2); 1904 if (plex->interpolated < 0) { 1905 ierr = DMPlexIsInterpolated_Internal(dm, &plex->interpolated);CHKERRQ(ierr); 1906 } else { 1907 #if defined (PETSC_USE_DEBUG) 1908 DMPlexInterpolatedFlag flg; 1909 1910 ierr = DMPlexIsInterpolated_Internal(dm, &flg);CHKERRQ(ierr); 1911 if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1912 #endif 1913 } 1914 *interpolated = plex->interpolated; 1915 PetscFunctionReturn(0); 1916 } 1917 1918 /*@ 1919 DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner). 1920 1921 Collective 1922 1923 Input Parameter: 1924 . dm - The DM object 1925 1926 Output Parameter: 1927 . interpolated - Flag whether the DM is interpolated 1928 1929 Level: intermediate 1930 1931 Notes: 1932 Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks. 1933 1934 This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks. 1935 1936 Developer Notes: 1937 Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID. 1938 1939 If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated. 1940 MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1941 if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED, 1942 otherwise sets plex->interpolatedCollective = plex->interpolated. 1943 1944 If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective. 1945 1946 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated() 1947 @*/ 1948 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1949 { 1950 DM_Plex *plex = (DM_Plex *) dm->data; 1951 PetscBool debug=PETSC_FALSE; 1952 PetscErrorCode ierr; 1953 1954 PetscFunctionBegin; 1955 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1956 PetscValidPointer(interpolated,2); 1957 ierr = PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);CHKERRQ(ierr); 1958 if (plex->interpolatedCollective < 0) { 1959 DMPlexInterpolatedFlag min, max; 1960 MPI_Comm comm; 1961 1962 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1963 ierr = DMPlexIsInterpolated(dm, &plex->interpolatedCollective);CHKERRQ(ierr); 1964 ierr = MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);CHKERRQ(ierr); 1965 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRQ(ierr); 1966 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1967 if (debug) { 1968 PetscMPIInt rank; 1969 1970 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1971 ierr = PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);CHKERRQ(ierr); 1972 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 1973 } 1974 } 1975 *interpolated = plex->interpolatedCollective; 1976 PetscFunctionReturn(0); 1977 } 1978