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 input 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 Developer Notes: 1259 It sets plex->interpolated = DMPLEX_INTERPOLATED_FULL. 1260 1261 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1262 @*/ 1263 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1264 { 1265 DMPlexInterpolatedFlag interpolated; 1266 DM idm, odm = dm; 1267 PetscSF sfPoint; 1268 PetscInt depth, dim, d; 1269 const char *name; 1270 PetscBool flg=PETSC_TRUE; 1271 PetscErrorCode ierr; 1272 1273 PetscFunctionBegin; 1274 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1275 PetscValidPointer(dmInt, 2); 1276 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1277 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1278 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1279 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1280 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1281 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1282 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1283 idm = dm; 1284 } else { 1285 for (d = 1; d < dim; ++d) { 1286 /* Create interpolated mesh */ 1287 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1288 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1289 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1290 if (depth > 0) { 1291 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1292 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1293 { 1294 PetscInt nroots; 1295 ierr = PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1296 if (nroots >= 0) {ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr);} 1297 } 1298 } 1299 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1300 odm = idm; 1301 } 1302 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1303 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1304 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1305 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr); 1306 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1307 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);} 1308 } 1309 { 1310 PetscBool isper; 1311 const PetscReal *maxCell, *L; 1312 const DMBoundaryType *bd; 1313 1314 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1315 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1316 } 1317 /* This function makes the mesh fully interpolated on all ranks */ 1318 { 1319 DM_Plex *plex = (DM_Plex *) idm->data; 1320 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1321 } 1322 *dmInt = idm; 1323 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1324 PetscFunctionReturn(0); 1325 } 1326 1327 /*@ 1328 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1329 1330 Collective on dmA 1331 1332 Input Parameter: 1333 . dmA - The DMPlex object with initial coordinates 1334 1335 Output Parameter: 1336 . dmB - The DMPlex object with copied coordinates 1337 1338 Level: intermediate 1339 1340 Note: This is typically used when adding pieces other than vertices to a mesh 1341 1342 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1343 @*/ 1344 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1345 { 1346 Vec coordinatesA, coordinatesB; 1347 VecType vtype; 1348 PetscSection coordSectionA, coordSectionB; 1349 PetscScalar *coordsA, *coordsB; 1350 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1351 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 1352 PetscBool lc = PETSC_FALSE; 1353 PetscErrorCode ierr; 1354 1355 PetscFunctionBegin; 1356 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1357 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1358 if (dmA == dmB) PetscFunctionReturn(0); 1359 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1360 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1361 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); 1362 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1363 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1364 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1365 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1366 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1367 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1368 if (!Nf) PetscFunctionReturn(0); 1369 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1370 if (!coordSectionB) { 1371 PetscInt dim; 1372 1373 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1374 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1375 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1376 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1377 } 1378 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1379 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1380 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1381 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1382 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1383 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); 1384 cS = cS - cStartA + cStartB; 1385 cE = vEndB; 1386 lc = PETSC_TRUE; 1387 } else { 1388 cS = vStartB; 1389 cE = vEndB; 1390 } 1391 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1392 for (v = vStartB; v < vEndB; ++v) { 1393 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1394 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1395 } 1396 if (lc) { /* localized coordinates */ 1397 PetscInt c; 1398 1399 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1400 PetscInt dof; 1401 1402 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1403 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1404 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1405 } 1406 } 1407 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1408 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1409 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1410 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1411 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1412 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1413 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1414 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1415 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1416 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1417 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1418 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1419 for (v = 0; v < vEndB-vStartB; ++v) { 1420 PetscInt offA, offB; 1421 1422 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1423 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1424 for (d = 0; d < spaceDim; ++d) { 1425 coordsB[offB+d] = coordsA[offA+d]; 1426 } 1427 } 1428 if (lc) { /* localized coordinates */ 1429 PetscInt c; 1430 1431 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1432 PetscInt dof, offA, offB; 1433 1434 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1435 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1436 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1437 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1438 } 1439 } 1440 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1441 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1442 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1443 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 /*@ 1448 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1449 1450 Collective on dm 1451 1452 Input Parameter: 1453 . dm - The complete DMPlex object 1454 1455 Output Parameter: 1456 . dmUnint - The DMPlex object with only cells and vertices 1457 1458 Level: intermediate 1459 1460 Notes: 1461 It does not copy over the coordinates. 1462 1463 Developer Notes: 1464 It sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1465 1466 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1467 @*/ 1468 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1469 { 1470 DMPlexInterpolatedFlag interpolated; 1471 DM udm; 1472 PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone; 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1477 PetscValidPointer(dmUnint, 2); 1478 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1479 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1480 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1481 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1482 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1483 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1484 *dmUnint = dm; 1485 PetscFunctionReturn(0); 1486 } 1487 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1488 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1489 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1490 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1491 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1492 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1493 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1494 for (c = cStart; c < cEnd; ++c) { 1495 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1496 1497 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1498 for (cl = 0; cl < closureSize*2; cl += 2) { 1499 const PetscInt p = closure[cl]; 1500 1501 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1502 } 1503 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1504 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1505 maxConeSize = PetscMax(maxConeSize, coneSize); 1506 } 1507 ierr = DMSetUp(udm);CHKERRQ(ierr); 1508 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1509 for (c = cStart; c < cEnd; ++c) { 1510 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1511 1512 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1513 for (cl = 0; cl < closureSize*2; cl += 2) { 1514 const PetscInt p = closure[cl]; 1515 1516 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1517 } 1518 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1519 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1520 } 1521 ierr = PetscFree(cone);CHKERRQ(ierr); 1522 ierr = DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1523 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1524 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1525 /* Reduce SF */ 1526 { 1527 PetscSF sfPoint, sfPointUn; 1528 const PetscSFNode *remotePoints; 1529 const PetscInt *localPoints; 1530 PetscSFNode *remotePointsUn; 1531 PetscInt *localPointsUn; 1532 PetscInt vEnd, numRoots, numLeaves, l; 1533 PetscInt numLeavesUn = 0, n = 0; 1534 PetscErrorCode ierr; 1535 1536 /* Get original SF information */ 1537 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1538 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1539 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1540 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1541 /* Allocate space for cells and vertices */ 1542 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1543 /* Fill in leaves */ 1544 if (vEnd >= 0) { 1545 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1546 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1547 for (l = 0; l < numLeaves; l++) { 1548 if (localPoints[l] < vEnd) { 1549 localPointsUn[n] = localPoints[l]; 1550 remotePointsUn[n].rank = remotePoints[l].rank; 1551 remotePointsUn[n].index = remotePoints[l].index; 1552 ++n; 1553 } 1554 } 1555 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1556 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1557 } 1558 } 1559 { 1560 PetscBool isper; 1561 const PetscReal *maxCell, *L; 1562 const DMBoundaryType *bd; 1563 1564 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1565 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1566 } 1567 /* This function makes the mesh fully uninterpolated on all ranks */ 1568 { 1569 DM_Plex *plex = (DM_Plex *) udm->data; 1570 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1571 } 1572 *dmUnint = udm; 1573 PetscFunctionReturn(0); 1574 } 1575 1576 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1577 { 1578 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1579 MPI_Comm comm; 1580 PetscErrorCode ierr; 1581 1582 PetscFunctionBegin; 1583 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1584 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1585 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1586 1587 if (depth == dim) { 1588 *interpolated = DMPLEX_INTERPOLATED_FULL; 1589 if (!dim) goto finish; 1590 1591 /* Check points at height = dim are vertices (have no cones) */ 1592 ierr = DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);CHKERRQ(ierr); 1593 for (p=pStart; p<pEnd; p++) { 1594 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1595 if (coneSize) { 1596 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1597 goto finish; 1598 } 1599 } 1600 1601 /* Check points at height < dim have cones */ 1602 for (h=0; h<dim; h++) { 1603 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 1604 for (p=pStart; p<pEnd; p++) { 1605 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1606 if (!coneSize) { 1607 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1608 goto finish; 1609 } 1610 } 1611 } 1612 } else if (depth == 1) { 1613 *interpolated = DMPLEX_INTERPOLATED_NONE; 1614 } else { 1615 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1616 } 1617 finish: 1618 PetscFunctionReturn(0); 1619 } 1620 1621 /*@ 1622 DMPlexIsInterpolated - Find out to what extent the DMPlex is topologically interpolated. 1623 1624 Not Collective 1625 1626 Input Parameter: 1627 . dm - The DM object 1628 1629 Output Parameter: 1630 . interpolated - Flag whether the DM is interpolated 1631 1632 Level: intermediate 1633 1634 Notes: 1635 Unlike DMPlexIsInterpolatedCollective(), this is NOT collective 1636 so the results can be different on different ranks in special cases. 1637 However, DMPlexInterpolate() guarantees the result is the same on all. 1638 1639 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED. 1640 1641 Developer Notes: 1642 Initially, plex->interpolated = DMPLEX_INTERPOLATED_INVALID. 1643 1644 If plex->interpolated == DMPLEX_INTERPOLATED_INVALID, DMPlexIsInterpolated_Internal() is called. 1645 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1646 DMPLEX_INTERPOLATED_NONE, DMPLEX_INTERPOLATED_PARTIAL or DMPLEX_INTERPOLATED_FULL. 1647 1648 If plex->interpolated != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolated. 1649 1650 DMPlexInterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_FULL, 1651 and DMPlexUninterpolate() sets plex->interpolated = DMPLEX_INTERPOLATED_NONE. 1652 1653 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective() 1654 @*/ 1655 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1656 { 1657 DM_Plex *plex = (DM_Plex *) dm->data; 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBegin; 1661 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1662 PetscValidPointer(interpolated,2); 1663 if (plex->interpolated < 0) { 1664 ierr = DMPlexIsInterpolated_Internal(dm, &plex->interpolated);CHKERRQ(ierr); 1665 } else { 1666 #if defined (PETSC_USE_DEBUG) 1667 DMPlexInterpolatedFlag flg; 1668 1669 ierr = DMPlexIsInterpolated_Internal(dm, &flg);CHKERRQ(ierr); 1670 if (flg != plex->interpolated) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1671 #endif 1672 } 1673 *interpolated = plex->interpolated; 1674 PetscFunctionReturn(0); 1675 } 1676 1677 /*@ 1678 DMPlexIsInterpolatedCollective - Find out to what extent the DMPlex is topologically interpolated (in collective manner). 1679 1680 Collective 1681 1682 Input Parameter: 1683 . dm - The DM object 1684 1685 Output Parameter: 1686 . interpolated - Flag whether the DM is interpolated 1687 1688 Level: intermediate 1689 1690 Notes: 1691 Unlike DMPlexIsInterpolated(), this is collective so the results are guaranteed to be the same on all ranks. 1692 1693 This function will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks. 1694 1695 Developer Notes: 1696 Initially, plex->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID. 1697 1698 If plex->interpolatedCollective == DMPLEX_INTERPOLATED_INVALID, this function calls DMPlexIsInterpolated() which sets plex->interpolated. 1699 MPI_Allreduce() is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1700 if plex->interpolated varies on different ranks, plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED, 1701 otherwise sets plex->interpolatedCollective = plex->interpolated. 1702 1703 If plex->interpolatedCollective != DMPLEX_INTERPOLATED_INVALID, this function just returns plex->interpolatedCollective. 1704 1705 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated() 1706 @*/ 1707 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1708 { 1709 DM_Plex *plex = (DM_Plex *) dm->data; 1710 PetscBool debug=PETSC_FALSE; 1711 PetscErrorCode ierr; 1712 1713 PetscFunctionBegin; 1714 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1715 PetscValidPointer(interpolated,2); 1716 ierr = PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);CHKERRQ(ierr); 1717 if (plex->interpolatedCollective < 0) { 1718 DMPlexInterpolatedFlag min, max; 1719 MPI_Comm comm; 1720 1721 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1722 ierr = DMPlexIsInterpolated(dm, &plex->interpolatedCollective);CHKERRQ(ierr); 1723 ierr = MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);CHKERRQ(ierr); 1724 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRQ(ierr); 1725 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1726 if (debug) { 1727 PetscMPIInt rank; 1728 1729 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1730 ierr = PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);CHKERRQ(ierr); 1731 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 1732 } 1733 } 1734 *interpolated = plex->interpolatedCollective; 1735 PetscFunctionReturn(0); 1736 } 1737