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, faceSizeAllT = 0, numCellFacesH = 0, faceT = 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 /* First get normal cell face size (we now allow hybrid cells to meet normal cells on either hybrid or normal faces */ 306 if (pStart[cellDepth] < pMax) {ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr);} 307 308 ierr = DMPlexGetConeSize(dm, pMax, &coneSizeH);CHKERRQ(ierr); 309 ierr = DMPlexGetCone(dm, pMax, &cone);CHKERRQ(ierr); 310 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);CHKERRQ(ierr); 311 if (faceSize < 0) { 312 PetscInt *sizes, minv, maxv; 313 314 /* count vertices of hybrid and non-hybrid faces */ 315 ierr = PetscCalloc1(numCellFacesH, &sizes);CHKERRQ(ierr); 316 for (cf = 0; cf < numCellFacesT; ++cf) { /* These are the non-hybrid faces */ 317 const PetscInt *cellFace = &cellFaces[-cf*faceSize]; 318 PetscInt f; 319 320 for (f = 0; f < -faceSize; ++f) sizes[cf] += (cellFace[f] >= 0 ? 1 : 0); 321 } 322 ierr = PetscSortInt(numCellFacesT, sizes);CHKERRQ(ierr); 323 minv = sizes[0]; 324 maxv = sizes[PetscMax(numCellFacesT-1, 0)]; 325 if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for non-hybrid face %D != %D", minv, maxv); 326 faceSizeAllT = minv; 327 ierr = PetscArrayzero(sizes, numCellFacesH);CHKERRQ(ierr); 328 for (cf = numCellFacesT; cf < numCellFacesH; ++cf) { /* These are the hybrid faces */ 329 const PetscInt *cellFace = &cellFaces[-cf*faceSize]; 330 PetscInt f; 331 332 for (f = 0; f < -faceSize; ++f) sizes[cf-numCellFacesT] += (cellFace[f] >= 0 ? 1 : 0); 333 } 334 ierr = PetscSortInt(numCellFacesH - numCellFacesT, sizes);CHKERRQ(ierr); 335 minv = sizes[0]; 336 maxv = sizes[PetscMax(numCellFacesH - numCellFacesT-1, 0)]; 337 ierr = PetscFree(sizes);CHKERRQ(ierr); 338 if (minv != maxv) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Different number of vertices for hybrid face %D != %D", minv, maxv); 339 faceSizeAllH = minv; 340 if (!faceSizeAll) faceSizeAll = faceSizeAllT; 341 } else { /* the size of the faces in hybrid cells is the same */ 342 faceSizeAll = faceSizeAllH = faceSizeAllT = faceSize; 343 } 344 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, cone, &numCellFacesH, &numCellFacesT, &faceSize, &cellFaces);CHKERRQ(ierr); 345 } else if (pEnd[cellDepth] > pStart[cellDepth]) { 346 ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr); 347 faceSizeAllH = faceSizeAllT = faceSizeAll; 348 } 349 if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 350 351 /* With hybrid grids, we first iterate on hybrid cells and start numbering the non-hybrid faces 352 Then, faces for non-hybrid cells are numbered. 353 This is to guarantee consistent orientations (all 0) of all the points in the cone of the hybrid cells */ 354 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 355 for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) { 356 PetscInt start, end; 357 358 start = outerloop == 0 ? pMax : pStart[cellDepth]; 359 end = outerloop == 0 ? pEnd[cellDepth] : pMax; 360 for (c = start; c < end; ++c) { 361 const PetscInt *cellFaces; 362 PetscInt numCellFaces, faceSize, faceSizeInc, faceSizeCheck, cf; 363 364 if (c < pMax) { 365 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 366 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 367 faceSizeCheck = faceSizeAll; 368 } else { /* Hybrid cell */ 369 const PetscInt *cone; 370 PetscInt numCellFacesN, coneSize; 371 372 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 373 if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH); 374 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 375 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 376 if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c); 377 faceSize = PetscMax(faceSize, -faceSize); 378 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize); 379 numCellFaces = numCellFacesN; /* process only non-hybrid faces */ 380 faceSizeCheck = faceSizeAllT; 381 } 382 faceSizeInc = faceSize; 383 for (cf = 0; cf < numCellFaces; ++cf) { 384 const PetscInt *cellFace = &cellFaces[cf*faceSizeInc]; 385 PetscInt faceSizeH = faceSize; 386 PetscHashIJKLKey key; 387 PetscHashIter iter; 388 PetscBool missing; 389 390 if (faceSizeInc == 2) { 391 key.i = PetscMin(cellFace[0], cellFace[1]); 392 key.j = PetscMax(cellFace[0], cellFace[1]); 393 key.k = PETSC_MAX_INT; 394 key.l = PETSC_MAX_INT; 395 } else { 396 key.i = cellFace[0]; 397 key.j = cellFace[1]; 398 key.k = cellFace[2]; 399 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 400 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 401 } 402 /* this check is redundant for non-hybrid meshes */ 403 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); 404 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 405 if (missing) { 406 ierr = PetscHashIJKLIterSet(faceTable, iter, face++);CHKERRQ(ierr); 407 if (c >= pMax) ++faceT; 408 } 409 } 410 if (c < pMax) { 411 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 412 } else { 413 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);CHKERRQ(ierr); 414 } 415 } 416 } 417 pEnd[faceDepth] = face; 418 419 /* Second pass for hybrid meshes: number hybrid faces */ 420 for (c = pMax; c < pEnd[cellDepth]; ++c) { 421 const PetscInt *cellFaces, *cone; 422 PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize; 423 424 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 425 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 426 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 427 if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH); 428 faceSize = PetscMax(faceSize, -faceSize); 429 for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */ 430 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 431 PetscHashIJKLKey key; 432 PetscHashIter iter; 433 PetscBool missing; 434 PetscInt faceSizeH = faceSize; 435 436 if (faceSize == 2) { 437 key.i = PetscMin(cellFace[0], cellFace[1]); 438 key.j = PetscMax(cellFace[0], cellFace[1]); 439 key.k = PETSC_MAX_INT; 440 key.l = PETSC_MAX_INT; 441 } else { 442 key.i = cellFace[0]; 443 key.j = cellFace[1]; 444 key.k = cellFace[2]; 445 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 446 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 447 } 448 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); 449 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 450 if (missing) {ierr = PetscHashIJKLIterSet(faceTable, iter, face++);CHKERRQ(ierr);} 451 } 452 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 453 } 454 faceH = face - pEnd[faceDepth]; 455 if (faceH) { 456 if (fMax == PETSC_DETERMINE) fMax = pEnd[faceDepth]; 457 else if (eMax == PETSC_DETERMINE) eMax = pEnd[faceDepth]; 458 else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of unassigned hybrid facets %D for cellDim %D and dimension %D", faceH, cellDim, dim); 459 } 460 pEnd[faceDepth] = face; 461 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 462 /* Count new points */ 463 for (d = 0; d <= depth; ++d) { 464 numPoints += pEnd[d]-pStart[d]; 465 } 466 ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr); 467 /* Set cone sizes */ 468 for (d = 0; d <= depth; ++d) { 469 PetscInt coneSize, p; 470 471 if (d == faceDepth) { 472 /* Now we have two cases: */ 473 if (faceSizeAll == faceSizeAllT) { 474 /* I see no way to do this if we admit faces of different shapes */ 475 for (p = pStart[d]; p < pEnd[d]-faceH; ++p) { 476 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 477 } 478 for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) { 479 ierr = DMPlexSetConeSize(idm, p, faceSizeAllH);CHKERRQ(ierr); 480 } 481 } else if (faceSizeAll == faceSizeAllH) { 482 for (p = pStart[d]; p < pStart[d]+faceT; ++p) { 483 ierr = DMPlexSetConeSize(idm, p, faceSizeAllT);CHKERRQ(ierr); 484 } 485 for (p = pStart[d]+faceT; p < pEnd[d]-faceH; ++p) { 486 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 487 } 488 for (p = pEnd[d]-faceH; p < pEnd[d]; ++p) { 489 ierr = DMPlexSetConeSize(idm, p, faceSizeAllH);CHKERRQ(ierr); 490 } 491 } else SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent faces sizes N: %D T: %D H: %D", faceSizeAll, faceSizeAllT, faceSizeAllH); 492 } else if (d == cellDepth) { 493 for (p = pStart[d]; p < pEnd[d]; ++p) { 494 /* Number of cell faces may be different from number of cell vertices*/ 495 if (p < pMax) { 496 ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr); 497 } else { 498 ierr = DMPlexGetFacesHybrid_Internal(dm, cellDim, p, &coneSize, NULL, NULL, NULL);CHKERRQ(ierr); 499 } 500 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 501 } 502 } else { 503 for (p = pStart[d]; p < pEnd[d]; ++p) { 504 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 505 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 506 } 507 } 508 } 509 ierr = DMSetUp(idm);CHKERRQ(ierr); 510 /* Get face cones from subsets of cell vertices */ 511 if (faceSizeAll > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 512 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 513 for (d = depth; d > cellDepth; --d) { 514 const PetscInt *cone; 515 PetscInt p; 516 517 for (p = pStart[d]; p < pEnd[d]; ++p) { 518 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 519 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr); 520 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr); 521 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr); 522 } 523 } 524 for (outerloop = 0, face = pStart[faceDepth]; outerloop < 2; outerloop++) { 525 PetscInt start, end; 526 527 start = outerloop == 0 ? pMax : pStart[cellDepth]; 528 end = outerloop == 0 ? pEnd[cellDepth] : pMax; 529 for (c = start; c < end; ++c) { 530 const PetscInt *cellFaces; 531 PetscInt numCellFaces, faceSize, faceSizeInc, cf; 532 533 if (c < pMax) { 534 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 535 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 536 } else { 537 const PetscInt *cone; 538 PetscInt numCellFacesN, coneSize; 539 540 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 541 if (coneSize != coneSizeH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid coneSize %D != %D", coneSize, coneSizeH); 542 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 543 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 544 if (numCellFaces != numCellFacesH) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected numCellFaces %D != %D for hybrid cell %D", numCellFaces, numCellFacesH, c); 545 faceSize = PetscMax(faceSize, -faceSize); 546 if (faceSize > 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSize); 547 numCellFaces = numCellFacesN; /* process only non-hybrid faces */ 548 } 549 faceSizeInc = faceSize; 550 for (cf = 0; cf < numCellFaces; ++cf) { 551 const PetscInt *cellFace = &cellFaces[cf*faceSizeInc]; 552 PetscHashIJKLKey key; 553 PetscHashIter iter; 554 PetscBool missing; 555 556 if (faceSizeInc == 2) { 557 key.i = PetscMin(cellFace[0], cellFace[1]); 558 key.j = PetscMax(cellFace[0], cellFace[1]); 559 key.k = PETSC_MAX_INT; 560 key.l = PETSC_MAX_INT; 561 } else { 562 key.i = cellFace[0]; 563 key.j = cellFace[1]; 564 key.k = cellFace[2]; 565 key.l = faceSizeInc > 3 ? (cellFace[3] < 0 ? faceSize = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 566 ierr = PetscSortInt(faceSizeInc, (PetscInt *) &key);CHKERRQ(ierr); 567 } 568 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 569 if (missing) { 570 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 571 ierr = PetscHashIJKLIterSet(faceTable, iter, face);CHKERRQ(ierr); 572 ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr); 573 } else { 574 const PetscInt *cone; 575 PetscInt coneSize, ornt, i, j, f; 576 577 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr); 578 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 579 /* Orient face: Do not allow reverse orientation at the first vertex */ 580 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 581 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 582 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); 583 /* - First find the initial vertex */ 584 for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break; 585 /* - Try forward comparison */ 586 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break; 587 if (j == faceSize) { 588 if ((faceSize == 2) && (i == 1)) ornt = -2; 589 else ornt = i; 590 } else { 591 /* - Try backward comparison */ 592 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break; 593 if (j == faceSize) { 594 if (i == 0) ornt = -faceSize; 595 else ornt = -i; 596 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c); 597 } 598 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 599 } 600 } 601 if (c < pMax) { 602 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 603 } else { 604 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSizeH, NULL, NULL, NULL, NULL, &cellFaces);CHKERRQ(ierr); 605 } 606 } 607 } 608 /* Second pass for hybrid meshes: orient hybrid faces */ 609 for (c = pMax; c < pEnd[cellDepth]; ++c) { 610 const PetscInt *cellFaces, *cone; 611 PetscInt numCellFaces, numCellFacesN, faceSize, cf, coneSize; 612 613 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 614 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 615 ierr = DMPlexGetRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 616 if (numCellFaces != numCellFacesH) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unexpected hybrid numCellFaces %D != %D", numCellFaces, numCellFacesH); 617 faceSize = PetscMax(faceSize, -faceSize); 618 for (cf = numCellFacesN; cf < numCellFaces; ++cf) { /* These are the hybrid faces */ 619 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 620 PetscHashIJKLKey key; 621 PetscHashIter iter; 622 PetscBool missing; 623 PetscInt faceSizeH = faceSize; 624 625 if (faceSize == 2) { 626 key.i = PetscMin(cellFace[0], cellFace[1]); 627 key.j = PetscMax(cellFace[0], cellFace[1]); 628 key.k = PETSC_MAX_INT; 629 key.l = PETSC_MAX_INT; 630 } else { 631 key.i = cellFace[0]; 632 key.j = cellFace[1]; 633 key.k = cellFace[2]; 634 key.l = faceSize > 3 ? (cellFace[3] < 0 ? faceSizeH = 3, PETSC_MAX_INT : cellFace[3]) : PETSC_MAX_INT; 635 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 636 } 637 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); 638 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 639 if (missing) { 640 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 641 ierr = PetscHashIJKLIterSet(faceTable, iter, face);CHKERRQ(ierr); 642 ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr); 643 } else { 644 PetscInt fv[4] = {0, 1, 2, 3}; 645 const PetscInt *cone; 646 PetscInt coneSize, ornt, i, j, f; 647 PetscBool q2h = PETSC_FALSE; 648 649 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr); 650 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 651 /* Orient face: Do not allow reverse orientation at the first vertex */ 652 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 653 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 654 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); 655 /* Hybrid faces are stored as tensor products of edges, so to compare them to normal faces, we have to flip */ 656 if (faceSize == 4 && c >= pMax && faceSizeAll != faceSizeAllT && f < pEnd[faceDepth] - faceH) {q2h = PETSC_TRUE; fv[2] = 3; fv[3] = 2;} 657 /* - First find the initial vertex */ 658 for (i = 0; i < faceSizeH; ++i) if (cellFace[fv[0]] == cone[i]) break; 659 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 */ 660 /* - Try forward comparison */ 661 for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+j)%faceSizeH]) break; 662 if (j == faceSizeH) { 663 if ((faceSizeH == 2) && (i == 1)) ornt = -2; 664 else ornt = i; 665 } else { 666 /* - Try backward comparison */ 667 for (j = 0; j < faceSizeH; ++j) if (cellFace[fv[j]] != cone[(i+faceSizeH-j)%faceSizeH]) break; 668 if (j == faceSizeH) { 669 if (i == 0) ornt = -faceSizeH; 670 else ornt = -i; 671 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine orientation of face %D in cell %D", f, c); 672 } 673 } else { 674 /* when matching hybrid faces in 3D, only few cases are possible. 675 Face traversal however can no longer follow the usual convention, this seems a serious issue to me */ 676 PetscInt tquad_map[4][4] = { {PETSC_MIN_INT, 0,PETSC_MIN_INT,PETSC_MIN_INT}, 677 { -1,PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT}, 678 {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT, 1}, 679 {PETSC_MIN_INT,PETSC_MIN_INT, -2,PETSC_MIN_INT} }; 680 PetscInt i2; 681 682 /* find the second vertex */ 683 for (i2 = 0; i2 < faceSizeH; ++i2) if (cellFace[fv[1]] == cone[i2]) break; 684 switch (faceSizeH) { 685 case 2: 686 ornt = i ? -2 : 0; 687 break; 688 case 4: 689 ornt = tquad_map[i][i2]; 690 break; 691 default: 692 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled face size %D for face %D in cell %D", faceSizeH, f, c); 693 694 } 695 } 696 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 697 } 698 } 699 ierr = DMPlexRestoreRawFacesHybrid_Internal(dm, cellDim, coneSize, cone, &numCellFaces, &numCellFacesN, &faceSize, &cellFaces);CHKERRQ(ierr); 700 } 701 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]); 702 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 703 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 704 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 705 ierr = DMPlexSetHybridBounds(idm, cMax, fMax, eMax, vMax);CHKERRQ(ierr); 706 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 707 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 708 PetscFunctionReturn(0); 709 } 710 711 PetscErrorCode DMPlexOrientCell(DM dm, PetscInt p, PetscInt masterConeSize, const PetscInt masterCone[]) 712 { 713 PetscInt coneSize; 714 PetscInt start1=0; 715 PetscBool reverse1=PETSC_FALSE; 716 const PetscInt *cone=NULL; 717 PetscErrorCode ierr; 718 719 PetscFunctionBegin; 720 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 721 if (!coneSize) PetscFunctionReturn(0); /* do nothing for points with no cone */ 722 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 723 ierr = DMPlexFixFaceOrientations_Orient_Private(coneSize, masterConeSize, masterCone, cone, &start1, &reverse1);CHKERRQ(ierr); 724 #if defined(PETSC_USE_DEBUG) 725 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]); 726 #endif 727 ierr = DMPlexOrientCell_Internal(dm, p, start1, reverse1);CHKERRQ(ierr); 728 #if defined(PETSC_USE_DEBUG) 729 { 730 PetscInt c; 731 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 732 for (c = 0; c < 2; c++) { 733 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); 734 } 735 } 736 #endif 737 PetscFunctionReturn(0); 738 } 739 740 PetscErrorCode DMPlexOrientCell_Internal(DM dm, PetscInt p, PetscInt start1, PetscBool reverse1) 741 { 742 PetscInt i, j, k, maxConeSize, coneSize, coneConeSize, supportSize, supportConeSize; 743 PetscInt start0, start; 744 PetscBool reverse0, reverse; 745 PetscInt newornt; 746 const PetscInt *cone=NULL, *support=NULL, *supportCone=NULL, *ornts=NULL; 747 PetscInt *newcone=NULL, *newornts=NULL; 748 PetscErrorCode ierr; 749 750 PetscFunctionBegin; 751 if (!start1 && !reverse1) PetscFunctionReturn(0); 752 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 753 if (!coneSize) PetscFunctionReturn(0); /* do nothing for points with no cone */ 754 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 755 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 756 /* permute p's cone and orientations */ 757 ierr = DMPlexGetConeOrientation(dm, p, &ornts);CHKERRQ(ierr); 758 ierr = DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newcone);CHKERRQ(ierr); 759 ierr = DMGetWorkArray(dm, maxConeSize, MPIU_INT, &newornts);CHKERRQ(ierr); 760 ierr = DMPlexFixFaceOrientations_Permute_Private(coneSize, cone, start1, reverse1, newcone);CHKERRQ(ierr); 761 ierr = DMPlexFixFaceOrientations_Permute_Private(coneSize, ornts, start1, reverse1, newornts);CHKERRQ(ierr); 762 /* if direction of p (face) is flipped, flip also p's cone points (edges) */ 763 if (reverse1) { 764 for (i=0; i<coneSize; i++) { 765 ierr = DMPlexGetConeSize(dm, cone[i], &coneConeSize);CHKERRQ(ierr); 766 ierr = DMPlexFixFaceOrientations_Translate_Private(newornts[i], &start0, &reverse0);CHKERRQ(ierr); 767 ierr = DMPlexFixFaceOrientations_Combine_Private(coneConeSize, start0, reverse0, 1, PETSC_FALSE, &start, &reverse);CHKERRQ(ierr); 768 ierr = DMPlexFixFaceOrientations_TranslateBack_Private(coneConeSize, start, reverse, &newornts[i]);CHKERRQ(ierr); 769 } 770 } 771 ierr = DMPlexSetConeOrientation(dm, p, newornts);CHKERRQ(ierr); 772 /* fix oriention of p within cones of p's support points */ 773 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 774 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 775 for (j=0; j<supportSize; j++) { 776 ierr = DMPlexGetCone(dm, support[j], &supportCone);CHKERRQ(ierr); 777 ierr = DMPlexGetConeSize(dm, support[j], &supportConeSize);CHKERRQ(ierr); 778 ierr = DMPlexGetConeOrientation(dm, support[j], &ornts);CHKERRQ(ierr); 779 for (k=0; k<supportConeSize; k++) { 780 if (supportCone[k] != p) continue; 781 ierr = DMPlexFixFaceOrientations_Translate_Private(ornts[k], &start0, &reverse0);CHKERRQ(ierr); 782 ierr = DMPlexFixFaceOrientations_Combine_Private(coneSize, start0, reverse0, start1, reverse1, &start, &reverse);CHKERRQ(ierr); 783 ierr = DMPlexFixFaceOrientations_TranslateBack_Private(coneSize, start, reverse, &newornt);CHKERRQ(ierr); 784 ierr = DMPlexInsertConeOrientation(dm, support[j], k, newornt);CHKERRQ(ierr); 785 } 786 } 787 /* rewrite cone */ 788 ierr = DMPlexSetCone(dm, p, newcone);CHKERRQ(ierr); 789 ierr = DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newcone);CHKERRQ(ierr); 790 ierr = DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &newornts);CHKERRQ(ierr); 791 PetscFunctionReturn(0); 792 } 793 794 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 795 { 796 PetscInt nleaves; 797 PetscInt nranks; 798 const PetscMPIInt *ranks=NULL; 799 const PetscInt *roffset=NULL, *rmine=NULL, *rremote=NULL; 800 PetscInt n, o, r; 801 PetscErrorCode ierr; 802 803 PetscFunctionBegin; 804 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 805 nleaves = roffset[nranks]; 806 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr); 807 for (r=0; r<nranks; r++) { 808 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 809 - to unify order with the other side */ 810 o = roffset[r]; 811 n = roffset[r+1] - o; 812 ierr = PetscArraycpy(&(*rmine1)[o], &rmine[o], n);CHKERRQ(ierr); 813 ierr = PetscArraycpy(&(*rremote1)[o], &rremote[o], n);CHKERRQ(ierr); 814 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr); 815 } 816 PetscFunctionReturn(0); 817 } 818 819 PetscErrorCode DMPlexOrientInterface(DM dm) 820 { 821 PetscSF sf=NULL; 822 PetscInt (*roots)[2], (*leaves)[2]; 823 PetscMPIInt (*rootsRanks)[2], (*leavesRanks)[2]; 824 const PetscInt *locals=NULL; 825 const PetscSFNode *remotes=NULL; 826 PetscInt nroots, nleaves, p, c; 827 PetscInt nranks, n, o, r; 828 const PetscMPIInt *ranks=NULL; 829 const PetscInt *roffset=NULL; 830 PetscInt *rmine1=NULL, *rremote1=NULL; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 831 const PetscInt *cone=NULL; 832 PetscInt coneSize, ind0; 833 MPI_Comm comm; 834 PetscMPIInt rank; 835 PetscInt debug = 0; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 840 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes);CHKERRQ(ierr); 841 if (nroots < 0) PetscFunctionReturn(0); 842 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 843 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr); 844 #if defined(PETSC_USE_DEBUG) 845 ierr = DMViewFromOptions(dm, NULL, "-before_fix_dm_view");CHKERRQ(ierr); 846 ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr); 847 #endif 848 ierr = SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr); 849 ierr = PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks);CHKERRQ(ierr); 850 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 851 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 852 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Roots\n");CHKERRQ(ierr);} 853 for (p = 0; p < nroots; ++p) { 854 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 855 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 856 /* Translate all points to root numbering */ 857 for (c = 0; c < 2; c++) { 858 if (coneSize > 1) { 859 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr); 860 if (ind0 < 0) { 861 roots[p][c] = cone[c]; 862 rootsRanks[p][c] = rank; 863 } else { 864 roots[p][c] = remotes[ind0].index; 865 rootsRanks[p][c] = remotes[ind0].rank; 866 } 867 } else { 868 roots[p][c] = -1; 869 rootsRanks[p][c] = -1; 870 } 871 } 872 } 873 if (debug) { 874 for (p = 0; p < nroots; ++p) { 875 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 876 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 877 if (coneSize > 1) { 878 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); 879 } 880 } 881 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 882 } 883 for (p = 0; p < nroots; ++p) { 884 for (c = 0; c < 2; c++) { 885 leaves[p][c] = -2; 886 leavesRanks[p][c] = -2; 887 } 888 } 889 ierr = PetscSFBcastBegin(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 890 ierr = PetscSFBcastBegin(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 891 ierr = PetscSFBcastEnd(sf, MPIU_2INT, roots, leaves);CHKERRQ(ierr); 892 ierr = PetscSFBcastEnd(sf, MPI_2INT, rootsRanks, leavesRanks);CHKERRQ(ierr); 893 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 894 if (debug && rank == 0) {ierr = PetscSynchronizedPrintf(comm, "Referred leaves\n");CHKERRQ(ierr);} 895 for (p = 0; p < nroots; ++p) { 896 if (leaves[p][0] < 0) continue; 897 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 898 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 899 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);} 900 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])) { 901 PetscInt masterCone[2]; 902 /* Translate these two cone points back to leave numbering */ 903 for (c = 0; c < 2; c++) { 904 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]); 905 /* Find index of rank leavesRanks[p][c] among remote ranks */ 906 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */ 907 ierr = PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r);CHKERRQ(ierr); 908 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]); 909 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */ 910 o = roffset[r]; 911 n = roffset[r+1] - o; 912 ierr = PetscFindInt(leaves[p][c], n, &rremote1[o], &ind0);CHKERRQ(ierr); 913 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]); 914 /* Get the corresponding local point */ 915 masterCone[c] = rmine1[o+ind0];CHKERRQ(ierr); 916 } 917 if (debug) {ierr = PetscSynchronizedPrintf(comm, "[%d] %D: masterCone=[%D %D]\n", rank, p, masterCone[0], masterCone[1]);CHKERRQ(ierr);} 918 /* Vaclav's note: Here we only compare first 2 points of the cone. Full cone size would lead to stronger self-checking. */ 919 ierr = DMPlexOrientCell(dm, p, 2, masterCone);CHKERRQ(ierr); 920 } 921 } 922 #if defined(PETSC_USE_DEBUG) 923 ierr = DMViewFromOptions(dm, NULL, "-after_fix_dm_view");CHKERRQ(ierr); 924 for (r = 0; r < nleaves; ++r) { 925 p = locals[r]; 926 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 927 if (!coneSize) continue; 928 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 929 for (c = 0; c < 2; c++) { 930 ierr = PetscFindInt(cone[c], nleaves, locals, &ind0);CHKERRQ(ierr); 931 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]); 932 if (leaves[p][c] != remotes[ind0].index || leavesRanks[p][c] != remotes[ind0].rank) { 933 if (leavesRanks[p][c] == rank) { 934 PetscInt ind1; 935 ierr = PetscFindInt(leaves[p][c], nleaves, locals, &ind1);CHKERRQ(ierr); 936 if (ind1 < 0) { 937 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]); 938 } else { 939 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); 940 } 941 } else { 942 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]); 943 } 944 } 945 } 946 } 947 #endif 948 if (debug) {ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);} 949 ierr = PetscFree4(roots, leaves, rootsRanks, leavesRanks);CHKERRQ(ierr); 950 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr); 951 PetscFunctionReturn(0); 952 } 953 954 static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[]) 955 { 956 PetscInt idx; 957 PetscMPIInt rank; 958 PetscBool flg; 959 PetscErrorCode ierr; 960 961 PetscFunctionBegin; 962 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr); 963 if (!flg) PetscFunctionReturn(0); 964 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 965 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr); 966 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]%s %D %s %D\n", rank, idxname, idx, valname, a[idx]);CHKERRQ(ierr);} 967 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 968 PetscFunctionReturn(0); 969 } 970 971 static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[]) 972 { 973 PetscInt idx; 974 PetscMPIInt rank; 975 PetscBool flg; 976 PetscErrorCode ierr; 977 978 PetscFunctionBegin; 979 ierr = PetscOptionsHasName(NULL, NULL, opt, &flg);CHKERRQ(ierr); 980 if (!flg) PetscFunctionReturn(0); 981 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 982 ierr = PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name);CHKERRQ(ierr); 983 if (idxname) { 984 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);} 985 } else { 986 for (idx = 0; idx < n; ++idx) {ierr = PetscSynchronizedPrintf(comm, "[%d]rank %D index %D\n", rank, a[idx].rank, a[idx].index);CHKERRQ(ierr);} 987 } 988 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr); 989 PetscFunctionReturn(0); 990 } 991 992 static PetscErrorCode DMPlexMapToLocalPoint(PetscHMapIJ roothash, const PetscInt localPoints[], PetscMPIInt rank, PetscSFNode remotePoint, PetscInt *localPoint) 993 { 994 PetscErrorCode ierr; 995 996 PetscFunctionBegin; 997 if (remotePoint.rank == rank) { 998 *localPoint = remotePoint.index; 999 } else { 1000 PetscHashIJKey key; 1001 PetscInt root; 1002 1003 key.i = remotePoint.index; 1004 key.j = remotePoint.rank; 1005 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 1006 if (root >= 0) { 1007 *localPoint = localPoints[root]; 1008 } else PetscFunctionReturn(1); 1009 } 1010 PetscFunctionReturn(0); 1011 } 1012 1013 /*@ 1014 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the PointSF in parallel, following local interpolation 1015 1016 Collective on dm 1017 1018 Input Parameters: 1019 + dm - The interpolated DM 1020 - pointSF - The initial SF without interpolated points 1021 1022 Output Parameter: 1023 . pointSF - The SF including interpolated points 1024 1025 Level: intermediate 1026 1027 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 1028 1029 .seealso: DMPlexInterpolate(), DMPlexUninterpolate() 1030 @*/ 1031 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 1032 { 1033 /* 1034 Okay, the algorithm is: 1035 - Take each point in the overlap (root) 1036 - Look at the neighboring points in the overlap (candidates) 1037 - Send these candidate points to neighbors 1038 - Neighbor checks for edge between root and candidate 1039 - If edge is found, it replaces candidate point with edge point 1040 - Send back the overwritten candidates (claims) 1041 - Original guy checks for edges, different from original candidate, and gets its own edge 1042 - This pair is put into SF 1043 1044 We need a new algorithm that tolerates groups larger than 2. 1045 - Take each point in the overlap (root) 1046 - Find all collections of points in the overlap which make faces (do early join) 1047 - Send collections as candidates (add size as first number) 1048 - Make sure to send collection to all owners of all overlap points in collection 1049 - Neighbor check for face in collections 1050 - If face is found, it replaces candidate point with face point 1051 - Send back the overwritten candidates (claims) 1052 - Original guy checks for faces, different from original candidate, and gets its own face 1053 - This pair is put into SF 1054 */ 1055 PetscHMapI leafhash; 1056 PetscHMapIJ roothash; 1057 const PetscInt *localPoints, *rootdegree; 1058 const PetscSFNode *remotePoints; 1059 PetscSFNode *candidates, *candidatesRemote, *claims; 1060 PetscSection candidateSection, candidateSectionRemote, claimSection; 1061 PetscInt numLeaves, l, numRoots, r, candidatesSize, candidatesRemoteSize; 1062 PetscMPIInt size, rank; 1063 PetscHashIJKey key; 1064 PetscBool debug = PETSC_FALSE; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 ierr = PetscOptionsHasName(NULL, ((PetscObject) dm)->prefix, "-dmplex_interp_debug", &debug);CHKERRQ(ierr); 1069 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 1070 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1071 ierr = PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1072 if (size < 2 || numRoots < 0) PetscFunctionReturn(0); 1073 ierr = DMPlexGetOverlap(dm, &r);CHKERRQ(ierr); 1074 if (r) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1075 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_interp_pre_view");CHKERRQ(ierr); 1076 ierr = PetscObjectViewFromOptions((PetscObject) pointSF, NULL, "-petscsf_interp_pre_view");CHKERRQ(ierr); 1077 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1078 /* Build hashes of points in the SF for efficient lookup */ 1079 ierr = PetscHMapICreate(&leafhash);CHKERRQ(ierr); 1080 ierr = PetscHMapIJCreate(&roothash);CHKERRQ(ierr); 1081 for (l = 0; l < numLeaves; ++l) { 1082 key.i = remotePoints[l].index; 1083 key.j = remotePoints[l].rank; 1084 ierr = PetscHMapISet(leafhash, localPoints[l], l);CHKERRQ(ierr); 1085 ierr = PetscHMapIJSet(roothash, key, l);CHKERRQ(ierr); 1086 } 1087 /* Compute root degree to identify shared points */ 1088 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 1089 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 1090 ierr = IntArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-interp_root_degree_view", "Root degree", "point", "degree", numRoots, rootdegree);CHKERRQ(ierr); 1091 /* Build a section / SFNode array of candidate points (face bd points) in the cone(support(leaf)), 1092 where each candidate is defined by a set of remote points (roots) for the other points that define the face. */ 1093 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);CHKERRQ(ierr); 1094 ierr = PetscSectionSetChart(candidateSection, 0, numRoots);CHKERRQ(ierr); 1095 { 1096 PetscHMapIJ facehash; 1097 1098 ierr = PetscHMapIJCreate(&facehash);CHKERRQ(ierr); 1099 for (l = 0; l < numLeaves; ++l) { 1100 const PetscInt localPoint = localPoints[l]; 1101 const PetscInt *support; 1102 PetscInt supportSize, s; 1103 1104 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);CHKERRQ(ierr);} 1105 ierr = DMPlexGetSupportSize(dm, localPoint, &supportSize);CHKERRQ(ierr); 1106 ierr = DMPlexGetSupport(dm, localPoint, &support);CHKERRQ(ierr); 1107 for (s = 0; s < supportSize; ++s) { 1108 const PetscInt face = support[s]; 1109 const PetscInt *cone; 1110 PetscInt coneSize, c, f, root; 1111 PetscBool isFace = PETSC_TRUE; 1112 1113 /* Only add face once */ 1114 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);CHKERRQ(ierr);} 1115 key.i = localPoint; 1116 key.j = face; 1117 ierr = PetscHMapIJGet(facehash, key, &f);CHKERRQ(ierr); 1118 if (f >= 0) continue; 1119 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1120 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 1121 /* If a cone point does not map to leaves on any proc, then do not put face in SF */ 1122 for (c = 0; c < coneSize; ++c) { 1123 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);CHKERRQ(ierr);} 1124 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1125 if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;} 1126 } 1127 if (isFace) { 1128 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found shared face %D\n", rank, face);CHKERRQ(ierr);} 1129 ierr = PetscHMapIJSet(facehash, key, l);CHKERRQ(ierr); 1130 ierr = PetscSectionAddDof(candidateSection, localPoint, coneSize);CHKERRQ(ierr); 1131 } 1132 } 1133 } 1134 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1135 ierr = PetscHMapIJClear(facehash);CHKERRQ(ierr); 1136 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 1137 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 1138 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 1139 for (l = 0; l < numLeaves; ++l) { 1140 const PetscInt localPoint = localPoints[l]; 1141 const PetscInt *support; 1142 PetscInt supportSize, s, offset, idx = 0; 1143 1144 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local point %D\n", rank, localPoint);CHKERRQ(ierr);} 1145 ierr = PetscSectionGetOffset(candidateSection, localPoint, &offset);CHKERRQ(ierr); 1146 ierr = DMPlexGetSupportSize(dm, localPoint, &supportSize);CHKERRQ(ierr); 1147 ierr = DMPlexGetSupport(dm, localPoint, &support);CHKERRQ(ierr); 1148 for (s = 0; s < supportSize; ++s) { 1149 const PetscInt face = support[s]; 1150 const PetscInt *cone; 1151 PetscInt coneSize, c, f, root; 1152 PetscBool isFace = PETSC_TRUE; 1153 1154 /* Only add face once */ 1155 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Support point %D\n", rank, face);CHKERRQ(ierr);} 1156 key.i = localPoint; 1157 key.j = face; 1158 ierr = PetscHMapIJGet(facehash, key, &f);CHKERRQ(ierr); 1159 if (f >= 0) continue; 1160 ierr = DMPlexGetConeSize(dm, face, &coneSize);CHKERRQ(ierr); 1161 ierr = DMPlexGetCone(dm, face, &cone);CHKERRQ(ierr); 1162 /* If a cone point does not map to leaves on any proc, then do not put face in SF */ 1163 for (c = 0; c < coneSize; ++c) { 1164 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Cone point %D\n", rank, cone[c]);CHKERRQ(ierr);} 1165 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1166 if (!rootdegree[cone[c]] && (root < 0)) {isFace = PETSC_FALSE; break;} 1167 } 1168 if (isFace) { 1169 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding shared face %D at idx %D\n", rank, face, idx);CHKERRQ(ierr);} 1170 ierr = PetscHMapIJSet(facehash, key, l);CHKERRQ(ierr); 1171 candidates[offset+idx].rank = -1; 1172 candidates[offset+idx++].index = coneSize-1; 1173 for (c = 0; c < coneSize; ++c) { 1174 if (cone[c] == localPoint) continue; 1175 if (rootdegree[cone[c]]) { 1176 candidates[offset+idx].rank = rank; 1177 candidates[offset+idx++].index = cone[c]; 1178 } else { 1179 ierr = PetscHMapIGet(leafhash, cone[c], &root);CHKERRQ(ierr); 1180 if (root < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot locate local point %D in SF", cone[c]); 1181 candidates[offset+idx++] = remotePoints[root]; 1182 } 1183 } 1184 } 1185 } 1186 } 1187 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1188 ierr = PetscHMapIJDestroy(&facehash);CHKERRQ(ierr); 1189 ierr = PetscObjectViewFromOptions((PetscObject) candidateSection, NULL, "-petscsection_interp_candidate_view");CHKERRQ(ierr); 1190 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates);CHKERRQ(ierr); 1191 } 1192 /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1193 /* Note that this section is indexed by offsets into leaves, not by point number */ 1194 { 1195 PetscSF sfMulti, sfInverse, sfCandidates; 1196 PetscInt *remoteOffsets; 1197 1198 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1199 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 1200 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);CHKERRQ(ierr); 1201 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);CHKERRQ(ierr); 1202 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);CHKERRQ(ierr); 1203 ierr = PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);CHKERRQ(ierr); 1204 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 1205 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1206 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 1207 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 1208 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 1209 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1210 1211 ierr = PetscObjectViewFromOptions((PetscObject) candidateSectionRemote, NULL, "-petscsection_interp_candidate_remote_view");CHKERRQ(ierr); 1212 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote);CHKERRQ(ierr); 1213 } 1214 /* */ 1215 { 1216 PetscInt idx; 1217 /* There is a section point for every leaf attached to a given root point */ 1218 for (r = 0, idx = 0; r < numRoots; ++r) { 1219 PetscInt deg; 1220 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1221 PetscInt offset, dof, d; 1222 1223 ierr = PetscSectionGetDof(candidateSectionRemote, idx, &dof);CHKERRQ(ierr); 1224 ierr = PetscSectionGetOffset(candidateSectionRemote, idx, &offset);CHKERRQ(ierr); 1225 for (d = 0; d < dof; ++d) { 1226 const PetscInt sizeInd = offset+d; 1227 const PetscInt numPoints = candidatesRemote[sizeInd].index; 1228 const PetscInt *join = NULL; 1229 PetscInt points[1024], p, joinSize; 1230 1231 points[0] = r; 1232 for (p = 0; p < numPoints; ++p) { 1233 ierr = DMPlexMapToLocalPoint(roothash, localPoints, rank, candidatesRemote[offset+(++d)], &points[p+1]); 1234 if (ierr) {d += numPoints-1 - p; break;} /* We got a point not in our overlap */ 1235 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking local candidate %D\n", rank, points[p+1]);CHKERRQ(ierr);} 1236 } 1237 if (ierr) continue; 1238 ierr = DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1239 if (joinSize == 1) { 1240 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Adding face %D at idx %D\n", rank, join[0], sizeInd);CHKERRQ(ierr);} 1241 candidatesRemote[sizeInd].rank = rank; 1242 candidatesRemote[sizeInd].index = join[0]; 1243 } 1244 ierr = DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1245 } 1246 } 1247 } 1248 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1249 } 1250 /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1251 { 1252 PetscSF sfMulti, sfClaims, sfPointNew; 1253 PetscSFNode *remotePointsNew; 1254 PetscHMapI claimshash; 1255 PetscInt *remoteOffsets, *localPointsNew; 1256 PetscInt claimsSize, pStart, pEnd, root, numLocalNew, p, d; 1257 1258 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 1259 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);CHKERRQ(ierr); 1260 ierr = PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);CHKERRQ(ierr); 1261 ierr = PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 1262 ierr = PetscSectionGetStorageSize(claimSection, &claimsSize);CHKERRQ(ierr); 1263 ierr = PetscMalloc1(claimsSize, &claims);CHKERRQ(ierr); 1264 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1265 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 1266 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 1267 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1268 ierr = PetscObjectViewFromOptions((PetscObject) claimSection, NULL, "-petscsection_interp_claim_view");CHKERRQ(ierr); 1269 ierr = SFNodeArrayViewFromOptions(PetscObjectComm((PetscObject) dm), "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims);CHKERRQ(ierr); 1270 /* Walk the original section of local supports and add an SF entry for each updated item */ 1271 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 1272 for (p = 0; p < numRoots; ++p) { 1273 PetscInt dof, offset; 1274 1275 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Checking root for claims %D\n", rank, p);CHKERRQ(ierr);} 1276 ierr = PetscSectionGetDof(candidateSection, p, &dof);CHKERRQ(ierr); 1277 ierr = PetscSectionGetOffset(candidateSection, p, &offset);CHKERRQ(ierr); 1278 for (d = 0; d < dof;) { 1279 if (claims[offset+d].rank >= 0) { 1280 const PetscInt faceInd = offset+d; 1281 const PetscInt numPoints = candidates[faceInd].index; 1282 const PetscInt *join = NULL; 1283 PetscInt joinSize, points[1024], c; 1284 1285 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);} 1286 points[0] = p; 1287 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[0]);CHKERRQ(ierr);} 1288 for (c = 0, ++d; c < numPoints; ++c, ++d) { 1289 key.i = candidates[offset+d].index; 1290 key.j = candidates[offset+d].rank; 1291 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 1292 points[c+1] = localPoints[root]; 1293 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] point %D\n", rank, points[c+1]);CHKERRQ(ierr);} 1294 } 1295 ierr = DMPlexGetJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1296 if (joinSize == 1) { 1297 if (debug) {ierr = PetscSynchronizedPrintf(PetscObjectComm((PetscObject) dm), "[%d] Found local face %D\n", rank, join[0]);CHKERRQ(ierr);} 1298 ierr = PetscHMapISet(claimshash, join[0], faceInd);CHKERRQ(ierr); 1299 } 1300 ierr = DMPlexRestoreJoin(dm, numPoints+1, points, &joinSize, &join);CHKERRQ(ierr); 1301 } else d += claims[offset+d].index+1; 1302 } 1303 } 1304 if (debug) {ierr = PetscSynchronizedFlush(PetscObjectComm((PetscObject) dm), NULL);CHKERRQ(ierr);} 1305 /* Create new pointSF from hashed claims */ 1306 ierr = PetscHMapIGetSize(claimshash, &numLocalNew);CHKERRQ(ierr); 1307 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1308 ierr = PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);CHKERRQ(ierr); 1309 ierr = PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);CHKERRQ(ierr); 1310 for (p = 0; p < numLeaves; ++p) { 1311 localPointsNew[p] = localPoints[p]; 1312 remotePointsNew[p].index = remotePoints[p].index; 1313 remotePointsNew[p].rank = remotePoints[p].rank; 1314 } 1315 p = numLeaves; 1316 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 1317 ierr = PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);CHKERRQ(ierr); 1318 for (p = numLeaves; p < numLeaves + numLocalNew; ++p) { 1319 PetscInt offset; 1320 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &offset);CHKERRQ(ierr); 1321 remotePointsNew[p] = claims[offset]; 1322 } 1323 ierr = PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);CHKERRQ(ierr); 1324 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1325 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 1326 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 1327 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 1328 } 1329 ierr = PetscHMapIDestroy(&leafhash);CHKERRQ(ierr); 1330 ierr = PetscHMapIJDestroy(&roothash);CHKERRQ(ierr); 1331 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 1332 ierr = PetscSectionDestroy(&candidateSectionRemote);CHKERRQ(ierr); 1333 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 1334 ierr = PetscFree(candidates);CHKERRQ(ierr); 1335 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 1336 ierr = PetscFree(claims);CHKERRQ(ierr); 1337 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 1338 PetscFunctionReturn(0); 1339 } 1340 1341 /*@ 1342 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1343 1344 Collective on dm 1345 1346 Input Parameters: 1347 + dm - The DMPlex object with only cells and vertices 1348 - dmInt - The interpolated DM 1349 1350 Output Parameter: 1351 . dmInt - The complete DMPlex object 1352 1353 Level: intermediate 1354 1355 Notes: 1356 It does not copy over the coordinates. 1357 1358 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1359 @*/ 1360 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1361 { 1362 DM idm, odm = dm; 1363 PetscSF sfPoint; 1364 PetscInt depth, dim, d; 1365 const char *name; 1366 PetscBool flg=PETSC_TRUE; 1367 PetscErrorCode ierr; 1368 1369 PetscFunctionBegin; 1370 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1371 PetscValidPointer(dmInt, 2); 1372 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1373 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1374 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1375 if ((depth == dim) || (dim <= 1)) { 1376 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1377 idm = dm; 1378 } else { 1379 for (d = 1; d < dim; ++d) { 1380 /* Create interpolated mesh */ 1381 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 1382 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 1383 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 1384 if (depth > 0) { 1385 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 1386 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 1387 ierr = DMPlexInterpolatePointSF(idm, sfPoint);CHKERRQ(ierr); 1388 } 1389 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 1390 odm = idm; 1391 } 1392 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 1393 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 1394 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 1395 ierr = DMCopyLabels(dm, idm);CHKERRQ(ierr); 1396 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL);CHKERRQ(ierr); 1397 if (flg) {ierr = DMPlexOrientInterface(idm);CHKERRQ(ierr);} 1398 } 1399 { 1400 PetscBool isper; 1401 const PetscReal *maxCell, *L; 1402 const DMBoundaryType *bd; 1403 1404 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1405 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 1406 } 1407 *dmInt = idm; 1408 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 1409 PetscFunctionReturn(0); 1410 } 1411 1412 /*@ 1413 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1414 1415 Collective on dmA 1416 1417 Input Parameter: 1418 . dmA - The DMPlex object with initial coordinates 1419 1420 Output Parameter: 1421 . dmB - The DMPlex object with copied coordinates 1422 1423 Level: intermediate 1424 1425 Note: This is typically used when adding pieces other than vertices to a mesh 1426 1427 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 1428 @*/ 1429 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1430 { 1431 Vec coordinatesA, coordinatesB; 1432 VecType vtype; 1433 PetscSection coordSectionA, coordSectionB; 1434 PetscScalar *coordsA, *coordsB; 1435 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1436 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 1437 PetscBool lc = PETSC_FALSE; 1438 PetscErrorCode ierr; 1439 1440 PetscFunctionBegin; 1441 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1442 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1443 if (dmA == dmB) PetscFunctionReturn(0); 1444 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 1445 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 1446 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); 1447 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 1448 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 1449 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 1450 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 1451 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 1452 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 1453 if (!Nf) PetscFunctionReturn(0); 1454 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 1455 if (!coordSectionB) { 1456 PetscInt dim; 1457 1458 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 1459 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 1460 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 1461 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 1462 } 1463 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1464 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 1465 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 1466 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 1467 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1468 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); 1469 cS = cS - cStartA + cStartB; 1470 cE = vEndB; 1471 lc = PETSC_TRUE; 1472 } else { 1473 cS = vStartB; 1474 cE = vEndB; 1475 } 1476 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 1477 for (v = vStartB; v < vEndB; ++v) { 1478 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 1479 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 1480 } 1481 if (lc) { /* localized coordinates */ 1482 PetscInt c; 1483 1484 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1485 PetscInt dof; 1486 1487 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1488 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 1489 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 1490 } 1491 } 1492 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1493 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 1494 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 1495 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1496 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1497 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 1498 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 1499 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 1500 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 1501 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 1502 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1503 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1504 for (v = 0; v < vEndB-vStartB; ++v) { 1505 PetscInt offA, offB; 1506 1507 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 1508 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 1509 for (d = 0; d < spaceDim; ++d) { 1510 coordsB[offB+d] = coordsA[offA+d]; 1511 } 1512 } 1513 if (lc) { /* localized coordinates */ 1514 PetscInt c; 1515 1516 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 1517 PetscInt dof, offA, offB; 1518 1519 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 1520 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 1521 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 1522 ierr = PetscArraycpy(coordsB + offB,coordsA + offA,dof);CHKERRQ(ierr); 1523 } 1524 } 1525 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 1526 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1527 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 1528 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 /*@ 1533 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1534 1535 Collective on dm 1536 1537 Input Parameter: 1538 . dm - The complete DMPlex object 1539 1540 Output Parameter: 1541 . dmUnint - The DMPlex object with only cells and vertices 1542 1543 Level: intermediate 1544 1545 Notes: 1546 It does not copy over the coordinates. 1547 1548 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 1549 @*/ 1550 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1551 { 1552 DM udm; 1553 PetscInt dim, vStart, vEnd, cStart, cEnd, cMax, c, maxConeSize = 0, *cone; 1554 PetscErrorCode ierr; 1555 1556 PetscFunctionBegin; 1557 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1558 PetscValidPointer(dmUnint, 2); 1559 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1560 if (dim <= 1) { 1561 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 1562 *dmUnint = dm; 1563 PetscFunctionReturn(0); 1564 } 1565 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1566 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1567 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1568 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 1569 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 1570 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 1571 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 1572 for (c = cStart; c < cEnd; ++c) { 1573 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1574 1575 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1576 for (cl = 0; cl < closureSize*2; cl += 2) { 1577 const PetscInt p = closure[cl]; 1578 1579 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1580 } 1581 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1582 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 1583 maxConeSize = PetscMax(maxConeSize, coneSize); 1584 } 1585 ierr = DMSetUp(udm);CHKERRQ(ierr); 1586 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 1587 for (c = cStart; c < cEnd; ++c) { 1588 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1589 1590 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1591 for (cl = 0; cl < closureSize*2; cl += 2) { 1592 const PetscInt p = closure[cl]; 1593 1594 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1595 } 1596 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1597 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 1598 } 1599 ierr = PetscFree(cone);CHKERRQ(ierr); 1600 ierr = DMPlexSetHybridBounds(udm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1601 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 1602 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 1603 /* Reduce SF */ 1604 { 1605 PetscSF sfPoint, sfPointUn; 1606 const PetscSFNode *remotePoints; 1607 const PetscInt *localPoints; 1608 PetscSFNode *remotePointsUn; 1609 PetscInt *localPointsUn; 1610 PetscInt vEnd, numRoots, numLeaves, l; 1611 PetscInt numLeavesUn = 0, n = 0; 1612 PetscErrorCode ierr; 1613 1614 /* Get original SF information */ 1615 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1616 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 1617 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 1618 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1619 /* Allocate space for cells and vertices */ 1620 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 1621 /* Fill in leaves */ 1622 if (vEnd >= 0) { 1623 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 1624 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 1625 for (l = 0; l < numLeaves; l++) { 1626 if (localPoints[l] < vEnd) { 1627 localPointsUn[n] = localPoints[l]; 1628 remotePointsUn[n].rank = remotePoints[l].rank; 1629 remotePointsUn[n].index = remotePoints[l].index; 1630 ++n; 1631 } 1632 } 1633 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 1634 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 1635 } 1636 } 1637 { 1638 PetscBool isper; 1639 const PetscReal *maxCell, *L; 1640 const DMBoundaryType *bd; 1641 1642 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 1643 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 1644 } 1645 1646 *dmUnint = udm; 1647 PetscFunctionReturn(0); 1648 } 1649