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