1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/hashmapi.h> 3 #include <petsc/private/hashmapij.h> 4 5 /* HashIJKL */ 6 7 #include <petsc/private/hashmap.h> 8 9 typedef struct _PetscHashIJKLKey { PetscInt i, j, k, l; } PetscHashIJKLKey; 10 11 #define PetscHashIJKLKeyHash(key) \ 12 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i),PetscHashInt((key).j)), \ 13 PetscHashCombine(PetscHashInt((key).k),PetscHashInt((key).l))) 14 15 #define PetscHashIJKLKeyEqual(k1,k2) \ 16 (((k1).i==(k2).i) ? ((k1).j==(k2).j) ? ((k1).k==(k2).k) ? ((k1).l==(k2).l) : 0 : 0 : 0) 17 18 PETSC_HASH_MAP(HashIJKL, PetscHashIJKLKey, PetscInt, PetscHashIJKLKeyHash, PetscHashIJKLKeyEqual, -1) 19 20 21 /* 22 DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell 23 This assumes that the mesh is not interpolated from the depth of point p to the vertices 24 */ 25 PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 26 { 27 const PetscInt *cone = NULL; 28 PetscInt coneSize; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 33 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 34 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 35 ierr = DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);CHKERRQ(ierr); 36 PetscFunctionReturn(0); 37 } 38 39 /* 40 DMPlexRestoreFaces_Internal - Restores the array 41 */ 42 PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 43 { 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 if (faces) { ierr = DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);CHKERRQ(ierr); } 48 PetscFunctionReturn(0); 49 } 50 51 /* 52 DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone 53 */ 54 PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 55 { 56 PetscInt *facesTmp; 57 PetscInt maxConeSize, maxSupportSize; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 62 if (faces && coneSize) PetscValidIntPointer(cone,4); 63 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 64 if (faces) {ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);CHKERRQ(ierr);} 65 switch (dim) { 66 case 1: 67 switch (coneSize) { 68 case 2: 69 if (faces) { 70 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 71 *faces = facesTmp; 72 } 73 if (numFaces) *numFaces = 2; 74 if (faceSize) *faceSize = 1; 75 break; 76 default: 77 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 78 } 79 break; 80 case 2: 81 switch (coneSize) { 82 case 3: 83 if (faces) { 84 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 85 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 86 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 87 *faces = facesTmp; 88 } 89 if (numFaces) *numFaces = 3; 90 if (faceSize) *faceSize = 2; 91 break; 92 case 4: 93 /* Vertices follow right hand rule */ 94 if (faces) { 95 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 96 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 97 facesTmp[4] = cone[2]; facesTmp[5] = cone[3]; 98 facesTmp[6] = cone[3]; facesTmp[7] = cone[0]; 99 *faces = facesTmp; 100 } 101 if (numFaces) *numFaces = 4; 102 if (faceSize) *faceSize = 2; 103 break; 104 default: 105 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 106 } 107 break; 108 case 3: 109 switch (coneSize) { 110 case 3: 111 if (faces) { 112 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 113 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 114 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 115 *faces = facesTmp; 116 } 117 if (numFaces) *numFaces = 3; 118 if (faceSize) *faceSize = 2; 119 break; 120 case 4: 121 /* Vertices of first face follow right hand rule and normal points away from last vertex */ 122 if (faces) { 123 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; 124 facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1]; 125 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3]; 126 facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3]; 127 *faces = facesTmp; 128 } 129 if (numFaces) *numFaces = 4; 130 if (faceSize) *faceSize = 3; 131 break; 132 case 8: 133 /* 7--------6 134 /| /| 135 / | / | 136 4--------5 | 137 | | | | 138 | | | | 139 | 1--------2 140 | / | / 141 |/ |/ 142 0--------3 143 */ 144 if (faces) { 145 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */ 146 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */ 147 facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */ 148 facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */ 149 facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */ 150 facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */ 151 *faces = facesTmp; 152 } 153 if (numFaces) *numFaces = 6; 154 if (faceSize) *faceSize = 4; 155 break; 156 default: 157 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 158 } 159 break; 160 default: 161 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 /* 167 DMPlexGetRawFacesHybrid_Internal - Gets groups of vertices that correspond to faces for the given cone using hybrid ordering (prisms) 168 */ 169 static PetscErrorCode DMPlexGetRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[]) 170 { 171 PetscInt *facesTmp; 172 PetscInt maxConeSize, maxSupportSize; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 177 if (faces && coneSize) PetscValidIntPointer(cone,4); 178 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 179 if (faces) {ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);CHKERRQ(ierr);} 180 switch (dim) { 181 case 1: 182 switch (coneSize) { 183 case 2: 184 if (faces) { 185 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 186 *faces = facesTmp; 187 } 188 if (numFaces) *numFaces = 2; 189 if (numFacesNotH) *numFacesNotH = 2; 190 if (faceSize) *faceSize = 1; 191 break; 192 default: 193 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 194 } 195 break; 196 case 2: 197 switch (coneSize) { 198 case 4: 199 if (faces) { 200 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 201 facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; 202 facesTmp[4] = cone[0]; facesTmp[5] = cone[2]; 203 facesTmp[6] = cone[1]; facesTmp[7] = cone[3]; 204 *faces = facesTmp; 205 } 206 if (numFaces) *numFaces = 4; 207 if (numFacesNotH) *numFacesNotH = 2; 208 if (faceSize) *faceSize = 2; 209 break; 210 default: 211 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 212 } 213 break; 214 case 3: 215 switch (coneSize) { 216 case 6: /* triangular prism */ 217 if (faces) { 218 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = -1; /* Bottom */ 219 facesTmp[4] = cone[3]; facesTmp[5] = cone[4]; facesTmp[6] = cone[5]; facesTmp[7] = -1; /* Top */ 220 facesTmp[8] = cone[0]; facesTmp[9] = cone[1]; facesTmp[10] = cone[3]; facesTmp[11] = cone[4]; /* Back left */ 221 facesTmp[12] = cone[1]; facesTmp[13] = cone[2]; facesTmp[14] = cone[4]; facesTmp[15] = cone[5]; /* Back right */ 222 facesTmp[16] = cone[2]; facesTmp[17] = cone[0]; facesTmp[18] = cone[5]; facesTmp[19] = cone[3]; /* Front */ 223 *faces = facesTmp; 224 } 225 if (numFaces) *numFaces = 5; 226 if (numFacesNotH) *numFacesNotH = 2; 227 if (faceSize) *faceSize = -4; 228 break; 229 default: 230 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 231 } 232 break; 233 default: 234 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 235 } 236 PetscFunctionReturn(0); 237 } 238 239 static PetscErrorCode DMPlexRestoreRawFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[]) 240 { 241 PetscErrorCode ierr; 242 243 PetscFunctionBegin; 244 if (faces) { ierr = DMRestoreWorkArray(dm, 0, MPIU_INT, (void *) faces);CHKERRQ(ierr); } 245 PetscFunctionReturn(0); 246 } 247 248 static PetscErrorCode DMPlexGetFacesHybrid_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *numFacesNotH, PetscInt *faceSize, const PetscInt *faces[]) 249 { 250 const PetscInt *cone = NULL; 251 PetscInt coneSize; 252 PetscErrorCode ierr; 253 254 PetscFunctionBegin; 255 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 256 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 257 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 258 ierr = DMPlexGetRawFacesHybrid_Internal(dm, dim, coneSize, cone, numFaces, numFacesNotH, faceSize, faces);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 /* This interpolates faces for cells at some stratum */ 263 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm) 264 { 265 DMLabel subpointMap; 266 PetscHashIJKL faceTable; 267 PetscInt *pStart, *pEnd; 268 PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d; 269 PetscInt coneSizeH = 0, faceSizeAllH = 0, numCellFacesH = 0, faceH, pMax = -1, dim, outerloop; 270 PetscInt cMax, fMax, eMax, vMax; 271 PetscErrorCode ierr; 272 273 PetscFunctionBegin; 274 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 275 /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */ 276 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 277 if (subpointMap) ++cellDim; 278 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 279 ++depth; 280 ++cellDepth; 281 cellDim -= depth - cellDepth; 282 ierr = PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);CHKERRQ(ierr); 283 for (d = depth-1; d >= faceDepth; --d) { 284 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);CHKERRQ(ierr); 285 } 286 ierr = DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);CHKERRQ(ierr); 287 pEnd[faceDepth] = pStart[faceDepth]; 288 for (d = faceDepth-1; d >= 0; --d) { 289 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 290 } 291 cMax = fMax = eMax = vMax = PETSC_DETERMINE; 292 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 293 if (cellDim == dim) { 294 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 295 pMax = cMax; 296 } else if (cellDim == dim -1) { 297 ierr = DMPlexGetHybridBounds(dm, &cMax, &fMax, NULL, NULL);CHKERRQ(ierr); 298 pMax = fMax; 299 } 300 pMax = pMax < 0 ? pEnd[cellDepth] : pMax; 301 if (pMax < pEnd[cellDepth]) { 302 const PetscInt *cellFaces, *cone; 303 PetscInt numCellFacesT, faceSize, cf; 304 305 ierr = DMPlexGetConeSize(dm, pMax, &coneSizeH);CHKERRQ(ierr); 306 ierr = DMPlexGetCone(dm, pMax, &cone);CHKERRQ(ierr); 307 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);CHKERRQ(ierr); 308 if (faceSize < 0) { 309 PetscInt *sizes, minv, maxv; 310 311 /* count vertices of hybrid and non-hybrid faces */ 312 ierr = PetscCalloc1(numCellFacesH, &sizes);CHKERRQ(ierr); 313 for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */ 314 const PetscInt *cellFace = &cellFaces[-cf*faceSize]; 315 PetscInt f; 316 317 for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0); 318 } 319 ierr = PetscSortInt(numCellFacesT, sizes);CHKERRQ(ierr); 320 minv = sizes[0]; 321 maxv = sizes[PetscMax(numCellFacesT-1, 0)]; 322 if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv); 323 faceSizeAll = minv; 324 ierr = PetscMemzero(sizes, numCellFacesH*sizeof(PetscInt));CHKERRQ(ierr); 325 for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */ 326 const PetscInt *cellFace = &cellFaces[-cf*faceSize]; 327 PetscInt f; 328 329 for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0); 330 } 331 ierr = PetscSortInt(numCellFacesH - numCellFacesT, sizes);CHKERRQ(ierr); 332 minv = sizes[0]; 333 maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)]; 334 ierr = PetscFree(sizes);CHKERRQ(ierr); 335 if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv); 336 faceSizeAllH = minv; 337 } else { /* the size of the faces in hybrid cells is the same */ 338 faceSizeAll = faceSizeAllH = faceSize; 339 } 340 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);CHKERRQ(ierr); 341 } else if (pEnd[cellDepth] > pStart[cellDepth]) { 342 ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr); 343 } 344 if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 345 346 /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid faces 347 Then, faces for non-hybrid cells are numbered. 348 This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */ 349 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 350 for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) { 351 PetscInt start, end; 352 353 start = outerloop == 0 ? pMax : pStart[cellDepth]; 354 end = outerloop == 0 ? pEnd[cellDepth] : pMax; 355 for (c = start; c < end; ++c) { 356 const PetscInt *cellFaces; 357 PetscInt numCellFaces, faceSize, faceSizeInc, cf; 358 359 if (c < pMax) { 360 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 361 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 362 } else { /* Hybrid cell */ 363 const PetscInt *cone; 364 PetscInt numCellFacesN, coneSize; 365 366 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 367 if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH); 368 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 369 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 370 if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c); 371 faceSize = PetscMax(faceSize, -faceSize); 372 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize); 373 numCellFaces = numCellFacesN; /* process only non-hybrid faces */ 374 } 375 faceSizeInc = faceSize; 376 for (cf = 0; cf < numCellFaces; ++cf) { 377 const PetscInt *cellFace = &cellFaces[cf*faceSizeInc]; 378 PetscInt faceSizeH = faceSize; 379 PetscHashIJKLKey key; 380 PetscHashIter iter; 381 PetscBool missing; 382 383 if (faceSizeInc == 2) { 384 key.i = PetscMin(cellFace[0], cellFace[1]); 385 key.j = PetscMax(cellFace[0], cellFace[1]); 386 key.k = PETSC_MAX_INT; 387 key.l = PETSC_MAX_INT; 388 } else { 389 key.i = cellFace[0]; 390 key.j = cellFace[1]; 391 key.k = cellFace[2]; 392 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 393 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 394 } 395 /* this check is redundant for non-hybrid meshes */ 396 if (faceSizeH != faceSizeAll) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected number of vertices for face %D of point %D -> %D != %D", cf, c, faceSizeH, faceSizeAll); 397 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 398 if (missing) {ierr = PetscHashIJKLIterSet(faceTable, iter, face++);CHKERRQ(ierr);} 399 } 400 if (c < pMax) { 401 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 402 } else { 403 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);CHKERRQ(ierr); 404 } 405 } 406 } 407 pEnd[faceDepth] = face; 408 409 /* Second pass for hybrid meshes: number hybrid faces */ 410 for (c = pMax; c < pEnd[cellDepth]; ++c) { 411 const PetscInt *cellFaces, *cone; 412 PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize; 413 414 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 415 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 416 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 417 if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH); 418 faceSize = PetscMax(faceSize, -faceSize); 419 for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */ 420 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 421 PetscHashIJKLKey key; 422 PetscHashIter iter; 423 PetscBool missing; 424 PetscInt faceSizeH = faceSize; 425 426 if (faceSize == 2) { 427 key.i = PetscMin(cellFace[0], cellFace[1]); 428 key.j = PetscMax(cellFace[0], cellFace[1]); 429 key.k = PETSC_MAX_INT; 430 key.l = PETSC_MAX_INT; 431 } else { 432 key.i = cellFace[0]; 433 key.j = cellFace[1]; 434 key.k = cellFace[2]; 435 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 436 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 437 } 438 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); 439 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 440 if (missing) {ierr = PetscHashIJKLIterSet(faceTable, iter, face++);CHKERRQ(ierr);} 441 } 442 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 443 } 444 faceH = face - pEnd[faceDepth]; 445 if (faceH) { 446 if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth]; 447 else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth]; 448 else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim); 449 } 450 pEnd[faceDepth] = face; 451 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 452 /* Count new points */ 453 for (d = 0; d <= depth; ++d) { 454 numPoints += pEnd[d]-pStart[d]; 455 } 456 ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr); 457 /* Set cone sizes */ 458 for (d = 0; d <= depth; ++d) { 459 PetscInt coneSize, p; 460 461 if (d == faceDepth) { 462 /* I see no way to do this if we admit faces of different shapes */ 463 for (p = pStart[d]; p < pEnd[d]-faceH; ++p) { 464 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 465 } 466 for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) { 467 ierr = DMPlexSetConeSize(idm, p, faceSizeAllH);CHKERRQ(ierr); 468 } 469 } else if (d == cellDepth) { 470 for (p = pStart[d]; p < pEnd[d]; ++p) { 471 /* Number of cell faces may be different from number of cell vertices*/ 472 if (p < pMax) { 473 ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr); 474 } else { 475 ierr = DMPlexGetFacesHybrid_Internal(dm, cellDim, p, &coneSize, NULL, NULL, NULL);CHKERRQ(ierr); 476 } 477 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 478 } 479 } else { 480 for (p = pStart[d]; p < pEnd[d]; ++p) { 481 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 482 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 483 } 484 } 485 } 486 ierr = DMSetUp(idm);CHKERRQ(ierr); 487 /* Get face cones from subsets of cell vertices */ 488 if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 489 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 490 for (d = depth; d > cellDepth; --d) { 491 const PetscInt *cone; 492 PetscInt p; 493 494 for (p = pStart[d]; p < pEnd[d]; ++p) { 495 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 496 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr); 497 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr); 498 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr); 499 } 500 } 501 for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) { 502 PetscInt start, end; 503 504 start = outerloop == 0 ? pMax : pStart[cellDepth]; 505 end = outerloop == 0 ? pEnd[cellDepth] : pMax; 506 for (c = start; c < end; ++c) { 507 const PetscInt *cellFaces; 508 PetscInt numCellFaces, faceSize, faceSizeInc, cf; 509 510 if (c < pMax) { 511 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 512 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 513 } else { 514 const PetscInt *cone; 515 PetscInt numCellFacesN, coneSize; 516 517 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 518 if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH); 519 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 520 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 521 if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c); 522 faceSize = PetscMax(faceSize, -faceSize); 523 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize); 524 numCellFaces = numCellFacesN; /* process only non-hybrid faces */ 525 } 526 faceSizeInc = faceSize; 527 for (cf = 0; cf < numCellFaces; ++cf) { 528 const PetscInt *cellFace = &cellFaces[cf*faceSizeInc]; 529 PetscHashIJKLKey key; 530 PetscHashIter iter; 531 PetscBool missing; 532 533 if (faceSizeInc == 2) { 534 key.i = PetscMin(cellFace[0], cellFace[1]); 535 key.j = PetscMax(cellFace[0], cellFace[1]); 536 key.k = PETSC_MAX_INT; 537 key.l = PETSC_MAX_INT; 538 } else { 539 key.i = cellFace[0]; 540 key.j = cellFace[1]; 541 key.k = cellFace[2]; 542 key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 543 ierr = PetscSortInt(faceSizeInc, (PetscInt *) &key);CHKERRQ(ierr); 544 } 545 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 546 if (missing) { 547 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 548 ierr = PetscHashIJKLIterSet(faceTable, iter, face);CHKERRQ(ierr); 549 ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr); 550 } else { 551 const PetscInt *cone; 552 PetscInt coneSize, ornt, i, j, f; 553 554 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr); 555 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 556 /* Orient face: Do not allow reverse orientation at the first vertex */ 557 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 558 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 559 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); 560 /* - First find the initial vertex */ 561 for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break; 562 /* - Try forward comparison */ 563 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break; 564 if (j == faceSize) { 565 if ((faceSize == 2) && (i == 1)) ornt = -2; 566 else ornt = i; 567 } else { 568 /* - Try backward comparison */ 569 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break; 570 if (j == faceSize) { 571 if (i == 0) ornt = -faceSize; 572 else ornt = -i; 573 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation"); 574 } 575 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 576 } 577 } 578 if (c < pMax) { 579 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 580 } else { 581 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);CHKERRQ(ierr); 582 } 583 } 584 } 585 /* Second pass for hybrid meshes: orient hybrid faces */ 586 for (c = pMax; c < pEnd[cellDepth]; ++c) { 587 const PetscInt *cellFaces, *cone; 588 PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize; 589 590 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 591 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 592 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 593 if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH); 594 faceSize = PetscMax(faceSize, -faceSize); 595 for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */ 596 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 597 PetscHashIJKLKey key; 598 PetscHashIter iter; 599 PetscBool missing; 600 PetscInt faceSizeH = faceSize; 601 602 if (faceSize == 2) { 603 key.i = PetscMin(cellFace[0], cellFace[1]); 604 key.j = PetscMax(cellFace[0], cellFace[1]); 605 key.k = PETSC_MAX_INT; 606 key.l = PETSC_MAX_INT; 607 } else { 608 key.i = cellFace[0]; 609 key.j = cellFace[1]; 610 key.k = cellFace[2]; 611 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 612 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 613 } 614 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); 615 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 616 if (missing) { 617 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 618 ierr = PetscHashIJKLIterSet(faceTable, iter, face);CHKERRQ(ierr); 619 ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr); 620 } else { 621 const PetscInt *cone; 622 PetscInt coneSize, ornt, i, j, f; 623 624 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr); 625 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 626 /* Orient face: Do not allow reverse orientation at the first vertex */ 627 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 628 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 629 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); 630 /* - First find the initial vertex */ 631 for (i = 0; i < faceSizeH; ++i) if (cellFace[0] == cone[i]) break; 632 /* - Try forward comparison */ 633 for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+j)%faceSizeH]) break; 634 if (j == faceSizeH) { 635 if ((faceSizeH == 2) && (i == 1)) ornt = -2; 636 else ornt = i; 637 } else { 638 /* - Try backward comparison */ 639 for (j = 0; j < faceSizeH; ++j) if (cellFace[j] != cone[(i+faceSizeH-j)%faceSizeH]) break; 640 if (j == faceSizeH) { 641 if (i == 0) ornt = -faceSizeH; 642 else ornt = -i; 643 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation"); 644 } 645 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 646 } 647 } 648 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 649 } 650 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]); 651 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 652 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 653 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 654 ierr = DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);CHKERRQ(ierr); 655 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 656 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 PetscErrorCode DMPlexOrientCell(DM dm, PetscInt p, PetscInt masterConeSize, const PetscInt masterCone[]) 661 { 662 PetscInt coneSize; 663 PetscInt start1=0; 664 PetscBool reverse1=PETSC_FALSE; 665 const PetscInt *cone=NULL; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 670 if (!coneSize) PetscFunctionReturn(0); /* do nothing for points with no cone */ 671 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 672 ierr = DMPlexFixFaceOrientations_Orient_Private(coneSize, masterConeSize, masterCone, cone, &start1, &reverse1);CHKERRQ(ierr); 673 #if defined(PETSC_USE_DEBUG) 674 if (PetscUnlikely(cone[start1] != masterCone[0])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[0]", start1, cone[start1], masterCone[0]); 675 #endif 676 ierr = DMPlexOrientCell_Internal(dm, p, start1, reverse1);CHKERRQ(ierr); 677 #if defined(PETSC_USE_DEBUG) 678 { 679 PetscInt c; 680 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 681 for (c = 0; c < 2; c++) { 682 if (PetscUnlikely(cone[c] != masterCone[c])) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "The algorithm above is wrong as cone[%d] = %d != %d = masterCone[%d]", c, cone[c], masterCone[c], c); 683 } 684 } 685 #endif 686 PetscFunctionReturn(0); 687 } 688 689 PetscErrorCode DMPlexOrientCell_Internal(DM dm, PetscInt p, PetscInt start1, PetscBool reverse1) 690 { 691 PetscInt i, j, k, maxConeSize, coneSize, coneConeSize, supportSize, supportConeSize; 692 PetscInt start0, start; 693 PetscBool reverse0, reverse; 694 PetscInt newornt; 695 const PetscInt *cone=NULL, *support=NULL, *supportCone=NULL, *ornts=NULL; 696 PetscInt *newcone=NULL, *newornts=NULL; 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 if (!start1 && !reverse1) PetscFunctionReturn(0); 701 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 702 if (!coneSize) PetscFunctionReturn(0); /* do nothing for points with no cone */ 703 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 704 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 705 /* permute p's cone and orientations */ 706 ierr = DMPlexGetConeOrientation(dm, p, &ornts);CHKERRQ(ierr); 707 ierr = DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newcone);CHKERRQ(ierr); 708 ierr = DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newornts);CHKERRQ(ierr); 709 ierr = DMPlexFixFaceOrientations_Permute_Private(coneSize, cone, start1, reverse1, newcone);CHKERRQ(ierr); 710 ierr = DMPlexFixFaceOrientations_Permute_Private(coneSize, ornts, start1, reverse1, newornts);CHKERRQ(ierr); 711 /* if direction of p (face) is flipped, flip also p's cone points (edges) */ 712 if (reverse1) { 713 for (i=0; i<coneSize; i++) { 714 ierr = DMPlexGetConeSize(dm, cone[i], &coneConeSize);CHKERRQ(ierr); 715 ierr = DMPlexFixFaceOrientations_Translate_Private(newornts[i], &start0, &reverse0);CHKERRQ(ierr); 716 ierr = DMPlexFixFaceOrientations_Combine_Private(coneConeSize, start0, reverse0, 1, PETSC_FALSE, &start, &reverse);CHKERRQ(ierr); 717 ierr = DMPlexFixFaceOrientations_TranslateBack_Private(coneConeSize, start, reverse, &newornts[i]);CHKERRQ(ierr); 718 } 719 } 720 ierr = DMPlexSetConeOrientation(dm, p, newornts);CHKERRQ(ierr); 721 /* fix oriention of p within cones of p's support points */ 722 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 723 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 724 for (j=0; j<supportSize; j++) { 725 ierr = DMPlexGetCone(dm, support[j], &supportCone);CHKERRQ(ierr); 726 ierr = DMPlexGetConeSize(dm, support[j], &supportConeSize);CHKERRQ(ierr); 727 ierr = DMPlexGetConeOrientation(dm, support[j], &ornts);CHKERRQ(ierr); 728 for (k=0; k<supportConeSize; k++) { 729 if (supportCone[k] != p) continue; 730 ierr = DMPlexFixFaceOrientations_Translate_Private(ornts[k], &start0, &reverse0);CHKERRQ(ierr); 731 ierr = DMPlexFixFaceOrientations_Combine_Private(coneSize, start0, reverse0, start1, reverse1, &start, &reverse);CHKERRQ(ierr); 732 ierr = DMPlexFixFaceOrientations_TranslateBack_Private(coneSize, start, reverse, &newornt);CHKERRQ(ierr); 733 ierr = DMPlexInsertConeOrientation(dm, support[j], k, newornt);CHKERRQ(ierr); 734 } 735 } 736 /* rewrite cone */ 737 ierr = DMPlexSetCone(dm, p, newcone);CHKERRQ(ierr); 738 ierr = DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newcone);CHKERRQ(ierr); 739 ierr = DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newornts);CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 744 { 745 PetscInt nleaves; 746 PetscInt nranks; 747 const PetscMPIInt *ranks=NULL; 748 const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL; 749 PetscInt n, o, r; 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 ierr = PetscSFGetRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 754 nleaves = roffset[nranks]; 755 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr); 756 for (r=0; r<nranks; r++) { 757 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 758 - to unify order with the other side */ 759 o = roffset[r]; 760 n = roffset[r+1] - o; 761 ierr = PetscMemcpy(&(*rmine1)[o], &rmine[o], n*sizeof(PetscInt));CHKERRQ(ierr); 762 ierr = PetscMemcpy(&(*rremote1)[o], &rremote[o], n*sizeof(PetscInt));CHKERRQ(ierr); 763 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr); 764 } 765 PetscFunctionReturn(0); 766 } 767 768 PetscErrorCode DMPlexOrientInterface(DM dm) 769 { 770 PetscSF sf=NULL; 771 PetscInt (*roots)[2], (*leaves)[2]; 772 PetscMPIInt (*rootsRanks)[2], (*leavesRanks)[2]; 773 const PetscInt *locals=NULL; 774 const PetscSFNode *remotes=NULL; 775 PetscInt nroots, nleaves, p, c; 776 PetscInt nranks, n, o, r; 777 const PetscMPIInt *ranks=NULL; 778 const PetscInt *roffset=NULL; 779 PetscInt *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 780 const PetscInt *cone=NULL; 781 PetscInt coneSize, ind0; 782 MPI_Comm comm; 783 PetscMPIInt rank; 784 PetscInt debug = 0; 785 PetscErrorCode ierr; 786 787 PetscFunctionBegin; 788 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 789 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);CHKERRQ(ierr); 790 if (nroots < 0) PetscFunctionReturn(0); 791 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 792 ierr = PetscSFGetRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr); 793 #if defined(PETSC_USE_DEBUG) 794 ierr = DMViewFromOptions(dm, NULL, "-before_fix_dm_view");CHKERRQ(ierr); 795 ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr); 796 #endif 797 ierr = SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr); 798 ierr = PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);CHKERRQ(ierr); 799 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 800 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 801 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Roots\n");CHKERRQ(ierr);} 802 for (p = 0; p < nroots; ++p) { 803 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 804 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 805 /* Translate all points to root numbering */ 806 for (c = 0; c < 2; c++) { 807 if (coneSize > 1) { 808 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr); 809 if (ind0 < 0) { 810 roots[p][c] = cone[c]; 811 rootsRanks[p][c] = rank; 812 } else { 813 roots[p][c] = remotes[ind0].index; 814 rootsRanks[p][c] = remotes[ind0].rank; 815 } 816 } else { 817 roots[p][c] = -1; 818 rootsRanks[p][c] = -1; 819 } 820 } 821 } 822 if (debug) { 823 for (p = 0; p < nroots; ++p) { 824 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 825 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 826 if (coneSize > 1) { 827 ierr = PetscSynchronizedPrintf(comm, "[%d] %D: cone=[%D %D] roots=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], roots[p][0], roots[p][1], rootsRanks[p][0], rootsRanks[p][1]);CHKERRQ(ierr); 828 } 829 } 830 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 831 } 832 for (p = 0; p < nroots; ++p) { 833 for (c = 0; c < 2; c++) { 834 leaves[p][c] = -2; 835 leavesRanks[p][c] = -2; 836 } 837 } 838 ierr = PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 839 ierr = PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 840 ierr = PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 841 ierr = PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 842 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 843 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Referred leaves\n");CHKERRQ(ierr);} 844 for (p = 0; p < nroots; ++p) { 845 if (leaves[p][0] < 0) continue; 846 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 847 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 848 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] %D: cone=[%D %D] leaves=[%D %D] roots=[%D %D] leavesRanks=[%D %D] rootsRanks=[%D %D]\n", rank, p, cone[0], cone[1], leaves[p][0], leaves[p][1], roots[p][0], roots[p][1], leavesRanks[p][0], leavesRanks[p][1], rootsRanks[p][0], rootsRanks[p][1]);CHKERRQ(ierr);} 849 if ((leaves[p][0] != roots[p][0]) || (leaves[p][1] != roots[p][1]) || (leavesRanks[p][0] != rootsRanks[p][0]) || (leavesRanks[p][0] != rootsRanks[p][0])) { 850 PetscInt masterCone[2]; 851 /* Translate these two cone points back to leave numbering */ 852 for (c = 0; c < 2; c++) { 853 if (leavesRanks[p][c] == rank) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - remote rank of point %D is the same rank",leavesRanks[p][c]); 854 /* Find index of rank leavesRanks[p][c] among remote ranks */ 855 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */ 856 ierr = PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);CHKERRQ(ierr); 857 if (PetscUnlikely(r < 0)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "this should never happen - rank %D not found among remote ranks",leavesRanks[p][c]); 858 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */ 859 o = roffset[r]; 860 n = roffset[r+1] - o; 861 ierr = PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);CHKERRQ(ierr); 862 if (PetscUnlikely(ind0 < 0)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No cone point of %D is connected to (%D, %D) - it seems there is missing connection in point SF!",p,ranks[r],leaves[p][c]); 863 /* Get the corresponding local point */ 864 masterCone[c] = rmine1[o+ind0];CHKERRQ(ierr); 865 } 866 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] %D: masterCone=[%D %D]\n", rank, p, masterCone[0], masterCone[1]);CHKERRQ(ierr);} 867 /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */ 868 ierr = DMPlexOrientCell(dm, p, 2, masterCone);CHKERRQ(ierr); 869 } 870 } 871 #if defined(PETSC_USE_DEBUG) 872 ierr = DMViewFromOptions(dm, NULL, "-after_fix_dm_view");CHKERRQ(ierr); 873 for (r = 0; r < nleaves; ++r) { 874 p = locals[r]; 875 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 876 if (!coneSize) continue; 877 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 878 for (c = 0; c < 2; c++) { 879 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr); 880 if (ind0 < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point SF contains %D but is missing its cone point cone[%D] = %D!", p, c, cone[c]); 881 if (leaves[p][c] != remotes[ind0].index || leavesRanks[p][c] != remotes[ind0].rank) { 882 if (leavesRanks[p][c] == rank) { 883 PetscInt ind1; 884 ierr = PetscFindInt(leaves[p][c], nleaves, locals, &ind1);CHKERRQ(ierr); 885 if (ind1 < 0) { 886 SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). The latter was not even found among the local SF points - it is probably broken!", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]); 887 } else { 888 SETERRQ9(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced %D --> (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leaves[p][c], remotes[ind1].rank, remotes[ind1].index); 889 } 890 } else { 891 SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D = locals[%d]: cone[%D]=%D --> (%D, %D) differs from the enforced (%D, %D). Is the algorithm above or the point SF broken?", p, r, c, cone[c], remotes[ind0].rank, remotes[ind0].index, leavesRanks[p][c], leaves[p][c]); 892 } 893 } 894 } 895 } 896 #endif 897 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 898 ierr = PetscFree4(roots, leaves, rootsRanks, leavesRanks);CHKERRQ(ierr); 899 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 /* This interpolates the PointSF in parallel following local interpolation */ 904 static PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF, PetscInt depth) 905 { 906 PetscMPIInt size, rank; 907 PetscInt p, c, d, dof, offset; 908 PetscInt numLeaves, numRoots, candidatesSize, candidatesRemoteSize; 909 const PetscInt *localPoints; 910 const PetscSFNode *remotePoints; 911 PetscSFNode *candidates, *candidatesRemote, *claims; 912 PetscSection candidateSection, candidateSectionRemote, claimSection; 913 PetscHMapI leafhash; 914 PetscHMapIJ roothash; 915 PetscHashIJKey key; 916 PetscErrorCode ierr; 917 918 PetscFunctionBegin; 919 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 920 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 921 ierr = PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 922 if (size < 2 || numRoots < 0) PetscFunctionReturn(0); 923 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 924 /* Build hashes of points in the SF for efficient lookup */ 925 ierr = PetscHMapICreate(&leafhash);CHKERRQ(ierr); 926 ierr = PetscHMapIJCreate(&roothash);CHKERRQ(ierr); 927 for (p = 0; p < numLeaves; ++p) { 928 ierr = PetscHMapISet(leafhash, localPoints[p], p);CHKERRQ(ierr); 929 key.i = remotePoints[p].index; 930 key.j = remotePoints[p].rank; 931 ierr = PetscHMapIJSet(roothash, key, p);CHKERRQ(ierr); 932 } 933 /* Build a section / SFNode array of candidate points in the single-level adjacency of leaves, 934 where each candidate is defined by the root entry for the other vertex that defines the edge. */ 935 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);CHKERRQ(ierr); 936 ierr = PetscSectionSetChart(candidateSection, 0, numRoots);CHKERRQ(ierr); 937 { 938 PetscInt leaf, root, idx, a, *adj = NULL; 939 for (p = 0; p < numLeaves; ++p) { 940 PetscInt adjSize = PETSC_DETERMINE; 941 ierr = DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);CHKERRQ(ierr); 942 for (a = 0; a < adjSize; ++a) { 943 ierr = PetscHMapIGet(leafhash, adj[a], &leaf);CHKERRQ(ierr); 944 if (leaf >= 0) {ierr = PetscSectionAddDof(candidateSection, localPoints[p], 1);CHKERRQ(ierr);} 945 } 946 } 947 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 948 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 949 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 950 for (p = 0; p < numLeaves; ++p) { 951 PetscInt adjSize = PETSC_DETERMINE; 952 ierr = PetscSectionGetOffset(candidateSection, localPoints[p], &offset);CHKERRQ(ierr); 953 ierr = DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);CHKERRQ(ierr); 954 for (idx = 0, a = 0; a < adjSize; ++a) { 955 ierr = PetscHMapIGet(leafhash, adj[a], &root);CHKERRQ(ierr); 956 if (root >= 0) candidates[offset+idx++] = remotePoints[root]; 957 } 958 } 959 ierr = PetscFree(adj);CHKERRQ(ierr); 960 } 961 /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 962 { 963 PetscSF sfMulti, sfInverse, sfCandidates; 964 PetscInt *remoteOffsets; 965 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 966 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 967 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);CHKERRQ(ierr); 968 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);CHKERRQ(ierr); 969 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);CHKERRQ(ierr); 970 ierr = PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);CHKERRQ(ierr); 971 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 972 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 973 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 974 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 975 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 976 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 977 } 978 /* Walk local roots and check for each remote candidate whether we know all required points, 979 either from owning it or having a root entry in the point SF. If we do we place a claim 980 by replacing the vertex number with our edge ID. */ 981 { 982 PetscInt idx, root, joinSize, vertices[2]; 983 const PetscInt *rootdegree, *join = NULL; 984 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 985 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 986 /* Loop remote edge connections and put in a claim if both vertices are known */ 987 for (idx = 0, p = 0; p < numRoots; ++p) { 988 for (d = 0; d < rootdegree[p]; ++d) { 989 ierr = PetscSectionGetDof(candidateSectionRemote, idx, &dof);CHKERRQ(ierr); 990 ierr = PetscSectionGetOffset(candidateSectionRemote, idx, &offset);CHKERRQ(ierr); 991 for (c = 0; c < dof; ++c) { 992 /* We own both vertices, so we claim the edge by replacing vertex with edge */ 993 if (candidatesRemote[offset+c].rank == rank) { 994 vertices[0] = p; vertices[1] = candidatesRemote[offset+c].index; 995 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 996 if (joinSize == 1) candidatesRemote[offset+c].index = join[0]; 997 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 998 continue; 999 } 1000 /* If we own one vertex and share a root with the other, we claim it */ 1001 key.i = candidatesRemote[offset+c].index; 1002 key.j = candidatesRemote[offset+c].rank; 1003 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 1004 if (root >= 0) { 1005 vertices[0] = p; vertices[1] = localPoints[root]; 1006 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 1007 if (joinSize == 1) { 1008 candidatesRemote[offset+c].index = join[0]; 1009 candidatesRemote[offset+c].rank = rank; 1010 } 1011 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 1012 } 1013 } 1014 idx++; 1015 } 1016 } 1017 } 1018 /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1019 { 1020 PetscSF sfMulti, sfClaims, sfPointNew; 1021 PetscHMapI claimshash; 1022 PetscInt size, pStart, pEnd, root, joinSize, numLocalNew; 1023 PetscInt *remoteOffsets, *localPointsNew, vertices[2]; 1024 const PetscInt *join = NULL; 1025 PetscSFNode *remotePointsNew; 1026 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1027 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);CHKERRQ(ierr); 1028 ierr = PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);CHKERRQ(ierr); 1029 ierr = PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1030 ierr = PetscSectionGetStorageSize(claimSection, &size);CHKERRQ(ierr); 1031 ierr = PetscMalloc1(size, &claims);CHKERRQ(ierr); 1032 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1033 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1034 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1035 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1036 /* Walk the original section of local supports and add an SF entry for each updated item */ 1037 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1038 for (p = 0; p < numRoots; ++p) { 1039 ierr = PetscSectionGetDof(candidateSection, p, &dof);CHKERRQ(ierr); 1040 ierr = PetscSectionGetOffset(candidateSection, p, &offset);CHKERRQ(ierr); 1041 for (d = 0; d < dof; ++d) { 1042 if (candidates[offset+d].index != claims[offset+d].index) { 1043 key.i = candidates[offset+d].index; 1044 key.j = candidates[offset+d].rank; 1045 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 1046 if (root >= 0) { 1047 vertices[0] = p; vertices[1] = localPoints[root]; 1048 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 1049 if (joinSize == 1) {ierr = PetscHMapISet(claimshash, join[0], offset+d);CHKERRQ(ierr);} 1050 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 1051 } 1052 } 1053 } 1054 } 1055 /* Create new pointSF from hashed claims */ 1056 ierr = PetscHMapIGetSize(claimshash, &numLocalNew);CHKERRQ(ierr); 1057 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1058 ierr = PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);CHKERRQ(ierr); 1059 ierr = PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);CHKERRQ(ierr); 1060 for (p = 0; p < numLeaves; ++p) { 1061 localPointsNew[p] = localPoints[p]; 1062 remotePointsNew[p].index = remotePoints[p].index; 1063 remotePointsNew[p].rank = remotePoints[p].rank; 1064 } 1065 p = numLeaves; 1066 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1067 ierr = PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);CHKERRQ(ierr); 1068 for (p = numLeaves; p < numLeaves + numLocalNew; ++p) { 1069 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &offset);CHKERRQ(ierr); 1070 remotePointsNew[p] = claims[offset]; 1071 } 1072 ierr = PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);CHKERRQ(ierr); 1073 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1074 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1075 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1076 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1077 } 1078 ierr = PetscHMapIDestroy(&leafhash);CHKERRQ(ierr); 1079 ierr = PetscHMapIJDestroy(&roothash);CHKERRQ(ierr); 1080 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1081 ierr = PetscSectionDestroy(&candidateSectionRemote);CHKERRQ(ierr); 1082 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1083 ierr = PetscFree(candidates);CHKERRQ(ierr); 1084 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1085 ierr = PetscFree(claims);CHKERRQ(ierr); 1086 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1087 PetscFunctionReturn(0); 1088 } 1089 1090 /*@C 1091 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1092 1093 Collective on DM 1094 1095 Input Parameters: 1096 + dm - The DMPlex object with only cells and vertices 1097 - dmInt - The interpolated DM 1098 1099 Output Parameter: 1100 . dmInt - The complete DMPlex object 1101 1102 Level: intermediate 1103 1104 Notes: 1105 It does not copy over the coordinates. 1106 1107 .keywords: mesh 1108 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1109 @*/ 1110 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1111 { 1112 DM idm, odm = dm; 1113 PetscSF sfPoint; 1114 PetscInt depth, dim, d; 1115 const char *name; 1116 PetscBool flg=PETSC_TRUE; 1117 PetscErrorCode ierr; 1118 1119 PetscFunctionBegin; 1120 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1121 PetscValidPointer(dmInt, 2); 1122 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1123 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1124 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1125 if ((depth == dim) || (dim <= 1)) { 1126 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1127 idm = dm; 1128 } else { 1129 for (d = 1; d < dim; ++d) { 1130 /* Create interpolated mesh */ 1131 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1132 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1133 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1134 if (depth > 0) { 1135 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1136 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1137 ierr = DMPlexInterpolatePointSF(idm, sfPoint, depth);CHKERRQ(ierr); 1138 } 1139 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1140 odm = idm; 1141 } 1142 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1143 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1144 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1145 ierr = DMCopyLabels(dm, idm);CHKERRQ(ierr); 1146 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1147 if (flg) {ierr = DMPlexOrientInterface(idm);CHKERRQ(ierr);} 1148 } 1149 { 1150 PetscBool isper; 1151 const PetscReal *maxCell, *L; 1152 const DMBoundaryType *bd; 1153 1154 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1155 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1156 } 1157 *dmInt = idm; 1158 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 /*@ 1163 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1164 1165 Collective on DM 1166 1167 Input Parameter: 1168 . dmA - The DMPlex object with initial coordinates 1169 1170 Output Parameter: 1171 . dmB - The DMPlex object with copied coordinates 1172 1173 Level: intermediate 1174 1175 Note: This is typically used when adding pieces other than vertices to a mesh 1176 1177 .keywords: mesh 1178 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1179 @*/ 1180 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1181 { 1182 Vec coordinatesA, coordinatesB; 1183 VecType vtype; 1184 PetscSection coordSectionA, coordSectionB; 1185 PetscScalar *coordsA, *coordsB; 1186 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1187 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 1188 PetscBool lc = PETSC_FALSE; 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1193 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1194 if (dmA == dmB) PetscFunctionReturn(0); 1195 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1196 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1197 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); 1198 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1199 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1200 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1201 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1202 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1203 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1204 if (!Nf) PetscFunctionReturn(0); 1205 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1206 if (!coordSectionB) { 1207 PetscInt dim; 1208 1209 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1210 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1211 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1212 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1213 } 1214 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1215 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1216 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1217 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1218 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1219 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); 1220 cS = cS - cStartA + cStartB; 1221 cE = vEndB; 1222 lc = PETSC_TRUE; 1223 } else { 1224 cS = vStartB; 1225 cE = vEndB; 1226 } 1227 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1228 for (v = vStartB; v < vEndB; ++v) { 1229 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1230 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1231 } 1232 if (lc) { /* localized coordinates */ 1233 PetscInt c; 1234 1235 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1236 PetscInt dof; 1237 1238 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1239 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1240 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1241 } 1242 } 1243 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1244 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1245 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1246 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1247 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1248 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1249 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1250 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1251 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1252 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1253 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1254 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1255 for (v = 0; v < vEndB-vStartB; ++v) { 1256 PetscInt offA, offB; 1257 1258 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1259 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1260 for (d = 0; d < spaceDim; ++d) { 1261 coordsB[offB+d] = coordsA[offA+d]; 1262 } 1263 } 1264 if (lc) { /* localized coordinates */ 1265 PetscInt c; 1266 1267 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1268 PetscInt dof, offA, offB; 1269 1270 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1271 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1272 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1273 ierr = PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));CHKERRQ(ierr); 1274 } 1275 } 1276 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1277 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1278 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1279 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1280 PetscFunctionReturn(0); 1281 } 1282 1283 /*@ 1284 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1285 1286 Collective on DM 1287 1288 Input Parameter: 1289 . dm - The complete DMPlex object 1290 1291 Output Parameter: 1292 . dmUnint - The DMPlex object with only cells and vertices 1293 1294 Level: intermediate 1295 1296 Notes: 1297 It does not copy over the coordinates. 1298 1299 .keywords: mesh 1300 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1301 @*/ 1302 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1303 { 1304 DM udm; 1305 PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone; 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1310 PetscValidPointer(dmUnint, 2); 1311 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1312 if (dim <= 1) { 1313 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1314 *dmUnint = dm; 1315 PetscFunctionReturn(0); 1316 } 1317 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1318 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1319 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1320 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1321 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1322 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1323 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1324 for (c = cStart; c < cEnd; ++c) { 1325 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1326 1327 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1328 for (cl = 0; cl < closureSize*2; cl += 2) { 1329 const PetscInt p = closure[cl]; 1330 1331 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1332 } 1333 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1334 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1335 maxConeSize = PetscMax(maxConeSize, coneSize); 1336 } 1337 ierr = DMSetUp(udm);CHKERRQ(ierr); 1338 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1339 for (c = cStart; c < cEnd; ++c) { 1340 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1341 1342 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1343 for (cl = 0; cl < closureSize*2; cl += 2) { 1344 const PetscInt p = closure[cl]; 1345 1346 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1347 } 1348 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1349 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1350 } 1351 ierr = PetscFree(cone);CHKERRQ(ierr); 1352 ierr = DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1353 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1354 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1355 /* Reduce SF */ 1356 { 1357 PetscSF sfPoint, sfPointUn; 1358 const PetscSFNode *remotePoints; 1359 const PetscInt *localPoints; 1360 PetscSFNode *remotePointsUn; 1361 PetscInt *localPointsUn; 1362 PetscInt vEnd, numRoots, numLeaves, l; 1363 PetscInt numLeavesUn = 0, n = 0; 1364 PetscErrorCode ierr; 1365 1366 /* Get original SF information */ 1367 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1368 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1369 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1370 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1371 /* Allocate space for cells and vertices */ 1372 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1373 /* Fill in leaves */ 1374 if (vEnd >= 0) { 1375 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1376 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1377 for (l = 0; l < numLeaves; l++) { 1378 if (localPoints[l] < vEnd) { 1379 localPointsUn[n] = localPoints[l]; 1380 remotePointsUn[n].rank = remotePoints[l].rank; 1381 remotePointsUn[n].index = remotePoints[l].index; 1382 ++n; 1383 } 1384 } 1385 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1386 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1387 } 1388 } 1389 { 1390 PetscBool isper; 1391 const PetscReal *maxCell, *L; 1392 const DMBoundaryType *bd; 1393 1394 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1395 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1396 } 1397 1398 *dmUnint = udm; 1399 PetscFunctionReturn(0); 1400 } 1401