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