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: intermediate 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 size, rank; 960 PetscHashIJKey key; 961 PetscBool debug = PETSC_FALSE; 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr); 966 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 967 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 968 ierr = PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 969 if (size < 2 || numRoots < 0) PetscFunctionReturn(0); 970 ierr = DMPlexGetOverlap(dm, &r);CHKERRQ(ierr); 971 if (r) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 972 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr); 973 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr); 974 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 975 /* Build hashes of points in the SF for efficient lookup */ 976 ierr = PetscHMapICreate(&leafhash);CHKERRQ(ierr); 977 ierr = PetscHMapIJCreate(&roothash);CHKERRQ(ierr); 978 for (l = 0; l < numLeaves; ++l) { 979 key.i = remotePoints[l].index; 980 key.j = remotePoints[l].rank; 981 ierr = PetscHMapISet(leafhash, localPoints[l], l);CHKERRQ(ierr); 982 ierr = PetscHMapIJSet(roothash, key, l);CHKERRQ(ierr); 983 } 984 /* Compute root degree to identify shared points */ 985 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 986 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 987 ierr = IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);CHKERRQ(ierr); 988 /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)), 989 where each candidate is defined by a set of remote points (roots) for the other points that define the face. */ 990 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);CHKERRQ(ierr); 991 ierr = PetscSectionSetChart(candidateSection, 0, numRoots);CHKERRQ(ierr); 992 { 993 PetscHMapIJ facehash; 994 995 ierr = PetscHMapIJCreate(&facehash);CHKERRQ(ierr); 996 for (l = 0; l < numLeaves; ++l) { 997 const PetscInt localPoint = localPoints[l]; 998 const PetscInt *support; 999 PetscInt supportSize, s; 1000 1001 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);CHKERRQ(ierr);} 1002 ierr = DMPlexGetSupportSize(dm, localPoint, &supportSize);CHKERRQ(ierr); 1003 ierr = DMPlexGetSupport(dm, localPoint, &support);CHKERRQ(ierr); 1004 for (s = 0; s < supportSize; ++s) { 1005 const PetscInt face = support[s]; 1006 const PetscInt *cone; 1007 PetscInt coneSize, c, f, root; 1008 PetscBool isFace = PETSC_TRUE; 1009 1010 /* Only add face once */ 1011 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);CHKERRQ(ierr);} 1012 key.i = localPoint; 1013 key.j = face; 1014 ierr = PetscHMapIJGet(facehash, key, &f);CHKERRQ(ierr); 1015 if (f >= 0) continue; 1016 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1017 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 1018 /* If a cone point does not map to leaves on any proc, then do not put face in SF */ 1019 for (c = 0; c < coneSize; ++c) { 1020 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);CHKERRQ(ierr);} 1021 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1022 if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;} 1023 } 1024 if (isFace) { 1025 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found shared face %D\n", rank, face);CHKERRQ(ierr);} 1026 ierr = PetscHMapIJSet(facehash, key, l);CHKERRQ(ierr); 1027 ierr = PetscSectionAddDof(candidateSection, localPoint, coneSize);CHKERRQ(ierr); 1028 } 1029 } 1030 } 1031 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1032 ierr = PetscHMapIJClear(facehash);CHKERRQ(ierr); 1033 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 1034 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 1035 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 1036 for (l = 0; l < numLeaves; ++l) { 1037 const PetscInt localPoint = localPoints[l]; 1038 const PetscInt *support; 1039 PetscInt supportSize, s, offset, idx = 0; 1040 1041 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);CHKERRQ(ierr);} 1042 ierr = PetscSectionGetOffset(candidateSection, localPoint, &offset);CHKERRQ(ierr); 1043 ierr = DMPlexGetSupportSize(dm, localPoint, &supportSize);CHKERRQ(ierr); 1044 ierr = DMPlexGetSupport(dm, localPoint, &support);CHKERRQ(ierr); 1045 for (s = 0; s < supportSize; ++s) { 1046 const PetscInt face = support[s]; 1047 const PetscInt *cone; 1048 PetscInt coneSize, c, f, root; 1049 PetscBool isFace = PETSC_TRUE; 1050 1051 /* Only add face once */ 1052 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);CHKERRQ(ierr);} 1053 key.i = localPoint; 1054 key.j = face; 1055 ierr = PetscHMapIJGet(facehash, key, &f);CHKERRQ(ierr); 1056 if (f >= 0) continue; 1057 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1058 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 1059 /* If a cone point does not map to leaves on any proc, then do not put face in SF */ 1060 for (c = 0; c < coneSize; ++c) { 1061 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);CHKERRQ(ierr);} 1062 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1063 if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;} 1064 } 1065 if (isFace) { 1066 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding shared face %D at idx %D\n", rank, face, idx);CHKERRQ(ierr);} 1067 ierr = PetscHMapIJSet(facehash, key, l);CHKERRQ(ierr); 1068 candidates[offset+idx].rank = -1; 1069 candidates[offset+idx++].index = coneSize-1; 1070 for (c = 0; c < coneSize; ++c) { 1071 if (cone[c] == localPoint) continue; 1072 if (rootdegree[cone[c]]) { 1073 candidates[offset+idx].rank = rank; 1074 candidates[offset+idx++].index = cone[c]; 1075 } else { 1076 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1077 if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]); 1078 candidates[offset+idx++] = remotePoints[root]; 1079 } 1080 } 1081 } 1082 } 1083 } 1084 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1085 ierr = PetscHMapIJDestroy(&facehash);CHKERRQ(ierr); 1086 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr); 1087 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr); 1088 } 1089 /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1090 /* Note that this section is indexed by offsets into leaves, not by point number */ 1091 { 1092 PetscSF sfMulti, sfInverse, sfCandidates; 1093 PetscInt *remoteOffsets; 1094 1095 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1096 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 1097 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);CHKERRQ(ierr); 1098 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);CHKERRQ(ierr); 1099 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);CHKERRQ(ierr); 1100 ierr = PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);CHKERRQ(ierr); 1101 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 1102 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1103 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1104 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 1105 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 1106 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1107 1108 ierr = PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr); 1109 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr); 1110 } 1111 /* */ 1112 { 1113 PetscInt idx; 1114 /* There is a section point for every leaf attached to a given root point */ 1115 for (r = 0, idx = 0; r < numRoots; ++r) { 1116 PetscInt deg; 1117 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1118 PetscInt offset, dof, d; 1119 1120 ierr = PetscSectionGetDof(candidateSectionRemote, idx, &dof);CHKERRQ(ierr); 1121 ierr = PetscSectionGetOffset(candidateSectionRemote, idx, &offset);CHKERRQ(ierr); 1122 for (d = 0; d < dof; ++d) { 1123 const PetscInt sizeInd = offset+d; 1124 const PetscInt numPoints = candidatesRemote[sizeInd].index; 1125 const PetscInt *join = NULL; 1126 PetscInt points[1024], p, joinSize; 1127 1128 points[0] = r; 1129 for (p = 0; p < numPoints; ++p) { 1130 ierr = DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]); 1131 if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */ 1132 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p+1]);CHKERRQ(ierr);} 1133 } 1134 if (ierr) continue; 1135 ierr = DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1136 if (joinSize == 1) { 1137 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding face %D at idx %D\n", rank, join[0], sizeInd);CHKERRQ(ierr);} 1138 candidatesRemote[sizeInd].rank = rank; 1139 candidatesRemote[sizeInd].index = join[0]; 1140 } 1141 ierr = DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1142 } 1143 } 1144 } 1145 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1146 } 1147 /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1148 { 1149 PetscSF sfMulti, sfClaims, sfPointNew; 1150 PetscSFNode *remotePointsNew; 1151 PetscHMapI claimshash; 1152 PetscInt *remoteOffsets, *localPointsNew; 1153 PetscInt claimsSize, pStart, pEnd, root, numLocalNew, p, d; 1154 1155 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1156 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);CHKERRQ(ierr); 1157 ierr = PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);CHKERRQ(ierr); 1158 ierr = PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1159 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr); 1160 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr); 1161 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1162 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1163 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1164 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1165 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr); 1166 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr); 1167 /* Walk the original section of local supports and add an SF entry for each updated item */ 1168 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1169 for (p = 0; p < numRoots; ++p) { 1170 PetscInt dof, offset; 1171 1172 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking root for claims %D\n", rank, p);CHKERRQ(ierr);} 1173 ierr = PetscSectionGetDof(candidateSection, p, &dof);CHKERRQ(ierr); 1174 ierr = PetscSectionGetOffset(candidateSection, p, &offset);CHKERRQ(ierr); 1175 for (d = 0; d < dof;) { 1176 if (claims[offset+d].rank >= 0) { 1177 const PetscInt faceInd = offset+d; 1178 const PetscInt numPoints = candidates[faceInd].index; 1179 const PetscInt *join = NULL; 1180 PetscInt joinSize, points[1024], c; 1181 1182 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);} 1183 points[0] = p; 1184 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);} 1185 for (c = 0, ++d; c < numPoints; ++c, ++d) { 1186 key.i = candidates[offset+d].index; 1187 key.j = candidates[offset+d].rank; 1188 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 1189 points[c+1] = localPoints[root]; 1190 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);} 1191 } 1192 ierr = DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1193 if (joinSize == 1) { 1194 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);} 1195 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr); 1196 } 1197 ierr = DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1198 } else d += claims[offset+d].index+1; 1199 } 1200 } 1201 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1202 /* Create new pointSF from hashed claims */ 1203 ierr = PetscHMapIGetSize(claimshash, &numLocalNew);CHKERRQ(ierr); 1204 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1205 ierr = PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);CHKERRQ(ierr); 1206 ierr = PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);CHKERRQ(ierr); 1207 for (p = 0; p < numLeaves; ++p) { 1208 localPointsNew[p] = localPoints[p]; 1209 remotePointsNew[p].index = remotePoints[p].index; 1210 remotePointsNew[p].rank = remotePoints[p].rank; 1211 } 1212 p = numLeaves; 1213 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1214 ierr = PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);CHKERRQ(ierr); 1215 for (p = numLeaves; p < numLeaves + numLocalNew; ++p) { 1216 PetscInt offset; 1217 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &offset);CHKERRQ(ierr); 1218 remotePointsNew[p] = claims[offset]; 1219 } 1220 ierr = PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);CHKERRQ(ierr); 1221 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1222 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1223 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1224 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1225 } 1226 ierr = PetscHMapIDestroy(&leafhash);CHKERRQ(ierr); 1227 ierr = PetscHMapIJDestroy(&roothash);CHKERRQ(ierr); 1228 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1229 ierr = PetscSectionDestroy(&candidateSectionRemote);CHKERRQ(ierr); 1230 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1231 ierr = PetscFree(candidates);CHKERRQ(ierr); 1232 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1233 ierr = PetscFree(claims);CHKERRQ(ierr); 1234 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 /*@ 1239 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1240 1241 Collective on dm 1242 1243 Input Parameters: 1244 + dm - The DMPlex object with only cells and vertices 1245 - dmInt - The interpolated DM 1246 1247 Output Parameter: 1248 . dmInt - The complete DMPlex object 1249 1250 Level: intermediate 1251 1252 Notes: 1253 It does not copy over the coordinates. 1254 1255 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1256 @*/ 1257 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1258 { 1259 DMPlexInterpolatedFlag interpolated; 1260 DM idm, odm = dm; 1261 PetscSF sfPoint; 1262 PetscInt depth, dim, d; 1263 const char *name; 1264 PetscBool flg=PETSC_TRUE; 1265 PetscErrorCode ierr; 1266 1267 PetscFunctionBegin; 1268 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1269 PetscValidPointer(dmInt, 2); 1270 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1271 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1272 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1273 if ((depth == dim) || (dim <= 1)) { 1274 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1275 idm = dm; 1276 } else { 1277 for (d = 1; d < dim; ++d) { 1278 /* Create interpolated mesh */ 1279 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1280 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1281 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1282 if (depth > 0) { 1283 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1284 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1285 ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr); 1286 } 1287 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1288 odm = idm; 1289 } 1290 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1291 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1292 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1293 ierr = DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE);CHKERRQ(ierr); 1294 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1295 if (flg) {ierr = DMPlexOrientInterface_Internal(idm);CHKERRQ(ierr);} 1296 } 1297 { 1298 PetscBool isper; 1299 const PetscReal *maxCell, *L; 1300 const DMBoundaryType *bd; 1301 1302 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1303 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1304 } 1305 *dmInt = idm; 1306 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 /*@ 1311 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1312 1313 Collective on dmA 1314 1315 Input Parameter: 1316 . dmA - The DMPlex object with initial coordinates 1317 1318 Output Parameter: 1319 . dmB - The DMPlex object with copied coordinates 1320 1321 Level: intermediate 1322 1323 Note: This is typically used when adding pieces other than vertices to a mesh 1324 1325 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1326 @*/ 1327 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1328 { 1329 Vec coordinatesA, coordinatesB; 1330 VecType vtype; 1331 PetscSection coordSectionA, coordSectionB; 1332 PetscScalar *coordsA, *coordsB; 1333 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1334 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 1335 PetscBool lc = PETSC_FALSE; 1336 PetscErrorCode ierr; 1337 1338 PetscFunctionBegin; 1339 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1340 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1341 if (dmA == dmB) PetscFunctionReturn(0); 1342 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1343 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1344 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); 1345 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1346 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1347 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1348 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1349 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1350 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1351 if (!Nf) PetscFunctionReturn(0); 1352 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1353 if (!coordSectionB) { 1354 PetscInt dim; 1355 1356 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1357 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1358 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1359 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1360 } 1361 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1362 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1363 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1364 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1365 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1366 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); 1367 cS = cS - cStartA + cStartB; 1368 cE = vEndB; 1369 lc = PETSC_TRUE; 1370 } else { 1371 cS = vStartB; 1372 cE = vEndB; 1373 } 1374 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1375 for (v = vStartB; v < vEndB; ++v) { 1376 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1377 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1378 } 1379 if (lc) { /* localized coordinates */ 1380 PetscInt c; 1381 1382 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1383 PetscInt dof; 1384 1385 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1386 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1387 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1388 } 1389 } 1390 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1391 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1392 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1393 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1394 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1395 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1396 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1397 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1398 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1399 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1400 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1401 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1402 for (v = 0; v < vEndB-vStartB; ++v) { 1403 PetscInt offA, offB; 1404 1405 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1406 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1407 for (d = 0; d < spaceDim; ++d) { 1408 coordsB[offB+d] = coordsA[offA+d]; 1409 } 1410 } 1411 if (lc) { /* localized coordinates */ 1412 PetscInt c; 1413 1414 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1415 PetscInt dof, offA, offB; 1416 1417 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1418 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1419 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1420 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1421 } 1422 } 1423 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1424 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1425 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1426 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1427 PetscFunctionReturn(0); 1428 } 1429 1430 /*@ 1431 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1432 1433 Collective on dm 1434 1435 Input Parameter: 1436 . dm - The complete DMPlex object 1437 1438 Output Parameter: 1439 . dmUnint - The DMPlex object with only cells and vertices 1440 1441 Level: intermediate 1442 1443 Notes: 1444 It does not copy over the coordinates. 1445 1446 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1447 @*/ 1448 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1449 { 1450 DM udm; 1451 PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone; 1452 PetscErrorCode ierr; 1453 1454 PetscFunctionBegin; 1455 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1456 PetscValidPointer(dmUnint, 2); 1457 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1458 if (dim <= 1) { 1459 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1460 *dmUnint = dm; 1461 PetscFunctionReturn(0); 1462 } 1463 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1464 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1465 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1466 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1467 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1468 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1469 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1470 for (c = cStart; c < cEnd; ++c) { 1471 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1472 1473 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1474 for (cl = 0; cl < closureSize*2; cl += 2) { 1475 const PetscInt p = closure[cl]; 1476 1477 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1478 } 1479 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1480 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1481 maxConeSize = PetscMax(maxConeSize, coneSize); 1482 } 1483 ierr = DMSetUp(udm);CHKERRQ(ierr); 1484 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1485 for (c = cStart; c < cEnd; ++c) { 1486 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1487 1488 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1489 for (cl = 0; cl < closureSize*2; cl += 2) { 1490 const PetscInt p = closure[cl]; 1491 1492 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1493 } 1494 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1495 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1496 } 1497 ierr = PetscFree(cone);CHKERRQ(ierr); 1498 ierr = DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1499 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1500 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1501 /* Reduce SF */ 1502 { 1503 PetscSF sfPoint, sfPointUn; 1504 const PetscSFNode *remotePoints; 1505 const PetscInt *localPoints; 1506 PetscSFNode *remotePointsUn; 1507 PetscInt *localPointsUn; 1508 PetscInt vEnd, numRoots, numLeaves, l; 1509 PetscInt numLeavesUn = 0, n = 0; 1510 PetscErrorCode ierr; 1511 1512 /* Get original SF information */ 1513 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1514 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1515 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1516 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1517 /* Allocate space for cells and vertices */ 1518 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1519 /* Fill in leaves */ 1520 if (vEnd >= 0) { 1521 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1522 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1523 for (l = 0; l < numLeaves; l++) { 1524 if (localPoints[l] < vEnd) { 1525 localPointsUn[n] = localPoints[l]; 1526 remotePointsUn[n].rank = remotePoints[l].rank; 1527 remotePointsUn[n].index = remotePoints[l].index; 1528 ++n; 1529 } 1530 } 1531 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1532 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1533 } 1534 } 1535 { 1536 PetscBool isper; 1537 const PetscReal *maxCell, *L; 1538 const DMBoundaryType *bd; 1539 1540 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1541 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1542 } 1543 1544 *dmUnint = udm; 1545 PetscFunctionReturn(0); 1546 } 1547 1548 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1549 { 1550 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1551 MPI_Comm comm; 1552 PetscErrorCode ierr; 1553 1554 PetscFunctionBegin; 1555 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1556 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1557 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1558 1559 if (depth == dim) { 1560 *interpolated = DMPLEX_INTERPOLATED_FULL; 1561 if (!dim) goto finish; 1562 1563 /* Check points at height = dim are vertices (have no cones) */ 1564 ierr = DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd);CHKERRQ(ierr); 1565 for (p=pStart; p<pEnd; p++) { 1566 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1567 if (coneSize) { 1568 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1569 goto finish; 1570 } 1571 } 1572 1573 /* Check points at height < dim have cones */ 1574 for (h=0; h<dim; h++) { 1575 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 1576 for (p=pStart; p<pEnd; p++) { 1577 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1578 if (!coneSize) { 1579 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1580 goto finish; 1581 } 1582 } 1583 } 1584 } else if (depth == 1) { 1585 *interpolated = DMPLEX_INTERPOLATED_NONE; 1586 } else { 1587 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1588 } 1589 finish: 1590 PetscFunctionReturn(0); 1591 } 1592 1593 /*@ 1594 DMPlexIsInterpolated - Find out whether this DM is interpolated, i.e. number of strata is equal to dimension. 1595 1596 Not Collective 1597 1598 Input Parameter: 1599 . dm - The DM object 1600 1601 Output Parameter: 1602 . interpolated - Flag whether the DM is interpolated 1603 1604 Level: intermediate 1605 1606 Notes: 1607 This is NOT collective so the results can be different on different ranks in special cases. 1608 However, DMPlexInterpolate() guarantees the result is the same on all. 1609 Unlike DMPlexIsInterpolatedCollective(), this cannot return DMPLEX_INTERPOLATED_MIXED. 1610 1611 .seealso: DMPlexInterpolate(), DMPlexIsInterpolatedCollective() 1612 @*/ 1613 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1614 { 1615 DM_Plex *plex = (DM_Plex *) dm->data; 1616 PetscErrorCode ierr; 1617 1618 PetscFunctionBegin; 1619 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1620 PetscValidPointer(interpolated,2); 1621 if (plex->interpolated < 0) { 1622 ierr = DMPlexIsInterpolated_Internal(dm, &plex->interpolated);CHKERRQ(ierr); 1623 } else { 1624 #if defined (PETSC_USE_DEBUG) 1625 DMPlexInterpolatedFlag flg; 1626 1627 ierr = DMPlexIsInterpolated_Internal(dm, &flg);CHKERRQ(ierr); 1628 if (flg != plex->interpolated) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "stashed DMPlexInterpolatedFlag is inconsistent"); 1629 #endif 1630 } 1631 *interpolated = plex->interpolated; 1632 PetscFunctionReturn(0); 1633 } 1634 1635 /*@ 1636 DMPlexIsInterpolatedCollective - Find out whether this DM is interpolated, i.e. number of strata is equal to dimension. 1637 1638 Not Collective 1639 1640 Input Parameter: 1641 . dm - The DM object 1642 1643 Output Parameter: 1644 . interpolated - Flag whether the DM is interpolated 1645 1646 Level: intermediate 1647 1648 Notes: 1649 This is collective so the results are always guaranteed to be the same on all ranks. 1650 Unlike DMPlexIsInterpolated(), this will return DMPLEX_INTERPOLATED_MIXED if the results of DMPlexIsInterpolated() are different on different ranks. 1651 1652 .seealso: DMPlexInterpolate(), DMPlexIsInterpolated() 1653 @*/ 1654 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1655 { 1656 DM_Plex *plex = (DM_Plex *) dm->data; 1657 PetscBool debug=PETSC_FALSE; 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBegin; 1661 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1662 PetscValidPointer(interpolated,2); 1663 ierr = PetscOptionsGetBool(((PetscObject) dm)->options, ((PetscObject) dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL);CHKERRQ(ierr); 1664 if (plex->interpolatedCollective < 0) { 1665 DMPlexInterpolatedFlag min, max; 1666 MPI_Comm comm; 1667 1668 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1669 ierr = DMPlexIsInterpolated(dm, &plex->interpolatedCollective);CHKERRQ(ierr); 1670 ierr = MPI_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm);CHKERRQ(ierr); 1671 ierr = MPI_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm);CHKERRQ(ierr); 1672 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1673 if (debug) { 1674 PetscMPIInt rank; 1675 1676 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1677 ierr = PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective]);CHKERRQ(ierr); 1678 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 1679 } 1680 if (PetscUnlikely(min < 0)) SETERRQ(comm, PETSC_ERR_PLIB, "Bug in DMPlexIsInterpolated: DMPlexInterpolatedFlag unitialized on some rank"); 1681 } 1682 *interpolated = plex->interpolatedCollective; 1683 PetscFunctionReturn(0); 1684 } 1685