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