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