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