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 */ 24 PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 25 { 26 const PetscInt *cone = NULL; 27 PetscInt maxConeSize, maxSupportSize, coneSize; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 32 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 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 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 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 63 if (faces) {ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), MPIU_INT, &facesTmp);CHKERRQ(ierr);} 64 switch (dim) { 65 case 1: 66 switch (coneSize) { 67 case 2: 68 if (faces) { 69 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 70 *faces = facesTmp; 71 } 72 if (numFaces) *numFaces = 2; 73 if (faceSize) *faceSize = 1; 74 break; 75 default: 76 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 77 } 78 break; 79 case 2: 80 switch (coneSize) { 81 case 3: 82 if (faces) { 83 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 84 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 85 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 86 *faces = facesTmp; 87 } 88 if (numFaces) *numFaces = 3; 89 if (faceSize) *faceSize = 2; 90 break; 91 case 4: 92 /* Vertices follow right hand rule */ 93 if (faces) { 94 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 95 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 96 facesTmp[4] = cone[2]; facesTmp[5] = cone[3]; 97 facesTmp[6] = cone[3]; facesTmp[7] = cone[0]; 98 *faces = facesTmp; 99 } 100 if (numFaces) *numFaces = 4; 101 if (faceSize) *faceSize = 2; 102 break; 103 default: 104 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 105 } 106 break; 107 case 3: 108 switch (coneSize) { 109 case 3: 110 if (faces) { 111 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 112 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 113 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 114 *faces = facesTmp; 115 } 116 if (numFaces) *numFaces = 3; 117 if (faceSize) *faceSize = 2; 118 break; 119 case 4: 120 /* Vertices of first face follow right hand rule and normal points away from last vertex */ 121 if (faces) { 122 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; 123 facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1]; 124 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3]; 125 facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3]; 126 *faces = facesTmp; 127 } 128 if (numFaces) *numFaces = 4; 129 if (faceSize) *faceSize = 3; 130 break; 131 case 8: 132 /* 7--------6 133 /| /| 134 / | / | 135 4--------5 | 136 | | | | 137 | | | | 138 | 1--------2 139 | / | / 140 |/ |/ 141 0--------3 142 */ 143 if (faces) { 144 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; /* Bottom */ 145 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; /* Top */ 146 facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */ 147 facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */ 148 facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */ 149 facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */ 150 *faces = facesTmp; 151 } 152 if (numFaces) *numFaces = 6; 153 if (faceSize) *faceSize = 4; 154 break; 155 default: 156 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 157 } 158 break; 159 default: 160 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 /* This interpolates faces for cells at some stratum */ 166 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm) 167 { 168 DMLabel subpointMap; 169 PetscHashIJKL faceTable; 170 PetscInt *pStart, *pEnd; 171 PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 176 /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */ 177 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 178 if (subpointMap) ++cellDim; 179 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 180 ++depth; 181 ++cellDepth; 182 cellDim -= depth - cellDepth; 183 ierr = PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);CHKERRQ(ierr); 184 for (d = depth-1; d >= faceDepth; --d) { 185 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);CHKERRQ(ierr); 186 } 187 ierr = DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);CHKERRQ(ierr); 188 pEnd[faceDepth] = pStart[faceDepth]; 189 for (d = faceDepth-1; d >= 0; --d) { 190 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 191 } 192 if (pEnd[cellDepth] > pStart[cellDepth]) {ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr);} 193 if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 194 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 195 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 196 const PetscInt *cellFaces; 197 PetscInt numCellFaces, faceSize, cf; 198 199 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 200 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 201 for (cf = 0; cf < numCellFaces; ++cf) { 202 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 203 PetscHashIJKLKey key; 204 PetscHashIter iter; 205 PetscBool missing; 206 207 if (faceSize == 2) { 208 key.i = PetscMin(cellFace[0], cellFace[1]); 209 key.j = PetscMax(cellFace[0], cellFace[1]); 210 key.k = 0; 211 key.l = 0; 212 } else { 213 key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0; 214 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 215 } 216 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 217 if (missing) {ierr = PetscHashIJKLIterSet(faceTable, iter, face++);CHKERRQ(ierr);} 218 } 219 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 220 } 221 pEnd[faceDepth] = face; 222 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 223 /* Count new points */ 224 for (d = 0; d <= depth; ++d) { 225 numPoints += pEnd[d]-pStart[d]; 226 } 227 ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr); 228 /* Set cone sizes */ 229 for (d = 0; d <= depth; ++d) { 230 PetscInt coneSize, p; 231 232 if (d == faceDepth) { 233 for (p = pStart[d]; p < pEnd[d]; ++p) { 234 /* I see no way to do this if we admit faces of different shapes */ 235 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 236 } 237 } else if (d == cellDepth) { 238 for (p = pStart[d]; p < pEnd[d]; ++p) { 239 /* Number of cell faces may be different from number of cell vertices*/ 240 ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr); 241 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 242 } 243 } else { 244 for (p = pStart[d]; p < pEnd[d]; ++p) { 245 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 246 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 247 } 248 } 249 } 250 ierr = DMSetUp(idm);CHKERRQ(ierr); 251 /* Get face cones from subsets of cell vertices */ 252 if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 253 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 254 for (d = depth; d > cellDepth; --d) { 255 const PetscInt *cone; 256 PetscInt p; 257 258 for (p = pStart[d]; p < pEnd[d]; ++p) { 259 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 260 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr); 261 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr); 262 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr); 263 } 264 } 265 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 266 const PetscInt *cellFaces; 267 PetscInt numCellFaces, faceSize, cf; 268 269 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 270 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 271 for (cf = 0; cf < numCellFaces; ++cf) { 272 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 273 PetscHashIJKLKey key; 274 PetscHashIter iter; 275 PetscBool missing; 276 277 if (faceSize == 2) { 278 key.i = PetscMin(cellFace[0], cellFace[1]); 279 key.j = PetscMax(cellFace[0], cellFace[1]); 280 key.k = 0; 281 key.l = 0; 282 } else { 283 key.i = cellFace[0]; 284 key.j = cellFace[1]; 285 key.k = cellFace[2]; 286 key.l = faceSize > 3 ? cellFace[3] : 0; 287 ierr = PetscSortInt(faceSize, (PetscInt *) &key);CHKERRQ(ierr); 288 } 289 ierr = PetscHashIJKLPut(faceTable, key, &iter, &missing);CHKERRQ(ierr); 290 if (missing) { 291 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 292 ierr = PetscHashIJKLIterSet(faceTable, iter, face);CHKERRQ(ierr); 293 ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr); 294 } else { 295 const PetscInt *cone; 296 PetscInt coneSize, ornt, i, j, f; 297 298 ierr = PetscHashIJKLIterGet(faceTable, iter, &f);CHKERRQ(ierr); 299 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 300 /* Orient face: Do not allow reverse orientation at the first vertex */ 301 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 302 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 303 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); 304 /* - First find the initial vertex */ 305 for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break; 306 /* - Try forward comparison */ 307 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break; 308 if (j == faceSize) { 309 if ((faceSize == 2) && (i == 1)) ornt = -2; 310 else ornt = i; 311 } else { 312 /* - Try backward comparison */ 313 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break; 314 if (j == faceSize) { 315 if (i == 0) ornt = -faceSize; 316 else ornt = -i; 317 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation"); 318 } 319 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 320 } 321 } 322 ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 323 } 324 if (face != pEnd[faceDepth]) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]); 325 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 326 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 327 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 328 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 329 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332 333 /* This interpolates the PointSF in parallel following local interpolation */ 334 static PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF, PetscInt depth) 335 { 336 PetscMPIInt size, rank; 337 PetscInt p, c, d, dof, offset; 338 PetscInt numLeaves, numRoots, candidatesSize, candidatesRemoteSize; 339 const PetscInt *localPoints; 340 const PetscSFNode *remotePoints; 341 PetscSFNode *candidates, *candidatesRemote, *claims; 342 PetscSection candidateSection, candidateSectionRemote, claimSection; 343 PetscHMapI leafhash; 344 PetscHMapIJ roothash; 345 PetscHashIJKey key; 346 PetscErrorCode ierr; 347 348 PetscFunctionBegin; 349 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 350 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 351 ierr = PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 352 if (size < 2 || numRoots < 0) PetscFunctionReturn(0); 353 ierr = PetscLogEventBegin(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 354 /* Build hashes of points in the SF for efficient lookup */ 355 ierr = PetscHMapICreate(&leafhash);CHKERRQ(ierr); 356 ierr = PetscHMapIJCreate(&roothash);CHKERRQ(ierr); 357 for (p = 0; p < numLeaves; ++p) { 358 ierr = PetscHMapISet(leafhash, localPoints[p], p);CHKERRQ(ierr); 359 key.i = remotePoints[p].index; 360 key.j = remotePoints[p].rank; 361 ierr = PetscHMapIJSet(roothash, key, p);CHKERRQ(ierr); 362 } 363 /* Build a section / SFNode array of candidate points in the single-level adjacency of leaves, 364 where each candidate is defined by the root entry for the other vertex that defines the edge. */ 365 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSection);CHKERRQ(ierr); 366 ierr = PetscSectionSetChart(candidateSection, 0, numRoots);CHKERRQ(ierr); 367 { 368 PetscInt leaf, root, idx, a, *adj = NULL; 369 for (p = 0; p < numLeaves; ++p) { 370 PetscInt adjSize = PETSC_DETERMINE; 371 ierr = DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);CHKERRQ(ierr); 372 for (a = 0; a < adjSize; ++a) { 373 ierr = PetscHMapIGet(leafhash, adj[a], &leaf);CHKERRQ(ierr); 374 if (leaf >= 0) {ierr = PetscSectionAddDof(candidateSection, localPoints[p], 1);CHKERRQ(ierr);} 375 } 376 } 377 ierr = PetscSectionSetUp(candidateSection);CHKERRQ(ierr); 378 ierr = PetscSectionGetStorageSize(candidateSection, &candidatesSize);CHKERRQ(ierr); 379 ierr = PetscMalloc1(candidatesSize, &candidates);CHKERRQ(ierr); 380 for (p = 0; p < numLeaves; ++p) { 381 PetscInt adjSize = PETSC_DETERMINE; 382 ierr = PetscSectionGetOffset(candidateSection, localPoints[p], &offset);CHKERRQ(ierr); 383 ierr = DMPlexGetAdjacency_Internal(dm, localPoints[p], PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &adjSize, &adj);CHKERRQ(ierr); 384 for (idx = 0, a = 0; a < adjSize; ++a) { 385 ierr = PetscHMapIGet(leafhash, adj[a], &root);CHKERRQ(ierr); 386 if (root >= 0) candidates[offset+idx++] = remotePoints[root]; 387 } 388 } 389 ierr = PetscFree(adj);CHKERRQ(ierr); 390 } 391 /* Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 392 { 393 PetscSF sfMulti, sfInverse, sfCandidates; 394 PetscInt *remoteOffsets; 395 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 396 ierr = PetscSFCreateInverseSF(sfMulti, &sfInverse);CHKERRQ(ierr); 397 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &candidateSectionRemote);CHKERRQ(ierr); 398 ierr = PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateSectionRemote);CHKERRQ(ierr); 399 ierr = PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateSectionRemote, &sfCandidates);CHKERRQ(ierr); 400 ierr = PetscSectionGetStorageSize(candidateSectionRemote, &candidatesRemoteSize);CHKERRQ(ierr); 401 ierr = PetscMalloc1(candidatesRemoteSize, &candidatesRemote);CHKERRQ(ierr); 402 ierr = PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 403 ierr = PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote);CHKERRQ(ierr); 404 ierr = PetscSFDestroy(&sfInverse);CHKERRQ(ierr); 405 ierr = PetscSFDestroy(&sfCandidates);CHKERRQ(ierr); 406 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 407 } 408 /* Walk local roots and check for each remote candidate whether we know all required points, 409 either from owning it or having a root entry in the point SF. If we do we place a claim 410 by replacing the vertex number with our edge ID. */ 411 { 412 PetscInt idx, root, joinSize, vertices[2]; 413 const PetscInt *rootdegree, *join = NULL; 414 ierr = PetscSFComputeDegreeBegin(pointSF, &rootdegree);CHKERRQ(ierr); 415 ierr = PetscSFComputeDegreeEnd(pointSF, &rootdegree);CHKERRQ(ierr); 416 /* Loop remote edge connections and put in a claim if both vertices are known */ 417 for (idx = 0, p = 0; p < numRoots; ++p) { 418 for (d = 0; d < rootdegree[p]; ++d) { 419 ierr = PetscSectionGetDof(candidateSectionRemote, idx, &dof);CHKERRQ(ierr); 420 ierr = PetscSectionGetOffset(candidateSectionRemote, idx, &offset);CHKERRQ(ierr); 421 for (c = 0; c < dof; ++c) { 422 /* We own both vertices, so we claim the edge by replacing vertex with edge */ 423 if (candidatesRemote[offset+c].rank == rank) { 424 vertices[0] = p; vertices[1] = candidatesRemote[offset+c].index; 425 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 426 if (joinSize == 1) candidatesRemote[offset+c].index = join[0]; 427 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 428 continue; 429 } 430 /* If we own one vertex and share a root with the other, we claim it */ 431 key.i = candidatesRemote[offset+c].index; 432 key.j = candidatesRemote[offset+c].rank; 433 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 434 if (root >= 0) { 435 vertices[0] = p; vertices[1] = localPoints[root]; 436 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 437 if (joinSize == 1) { 438 candidatesRemote[offset+c].index = join[0]; 439 candidatesRemote[offset+c].rank = rank; 440 } 441 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 442 } 443 } 444 idx++; 445 } 446 } 447 } 448 /* Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 449 { 450 PetscSF sfMulti, sfClaims, sfPointNew; 451 PetscHMapI claimshash; 452 PetscInt size, pStart, pEnd, root, joinSize, numLocalNew; 453 PetscInt *remoteOffsets, *localPointsNew, vertices[2]; 454 const PetscInt *join = NULL; 455 PetscSFNode *remotePointsNew; 456 ierr = PetscSFGetMultiSF(pointSF, &sfMulti);CHKERRQ(ierr); 457 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &claimSection);CHKERRQ(ierr); 458 ierr = PetscSFDistributeSection(sfMulti, candidateSectionRemote, &remoteOffsets, claimSection);CHKERRQ(ierr); 459 ierr = PetscSFCreateSectionSF(sfMulti, candidateSectionRemote, remoteOffsets, claimSection, &sfClaims);CHKERRQ(ierr); 460 ierr = PetscSectionGetStorageSize(claimSection, &size);CHKERRQ(ierr); 461 ierr = PetscMalloc1(size, &claims);CHKERRQ(ierr); 462 ierr = PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 463 ierr = PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims);CHKERRQ(ierr); 464 ierr = PetscSFDestroy(&sfClaims);CHKERRQ(ierr); 465 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 466 /* Walk the original section of local supports and add an SF entry for each updated item */ 467 ierr = PetscHMapICreate(&claimshash);CHKERRQ(ierr); 468 for (p = 0; p < numRoots; ++p) { 469 ierr = PetscSectionGetDof(candidateSection, p, &dof);CHKERRQ(ierr); 470 ierr = PetscSectionGetOffset(candidateSection, p, &offset);CHKERRQ(ierr); 471 for (d = 0; d < dof; ++d) { 472 if (candidates[offset+d].index != claims[offset+d].index) { 473 key.i = candidates[offset+d].index; 474 key.j = candidates[offset+d].rank; 475 ierr = PetscHMapIJGet(roothash, key, &root);CHKERRQ(ierr); 476 if (root >= 0) { 477 vertices[0] = p; vertices[1] = localPoints[root]; 478 ierr = DMPlexGetJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 479 if (joinSize == 1) {ierr = PetscHMapISet(claimshash, join[0], offset+d);CHKERRQ(ierr);} 480 ierr = DMPlexRestoreJoin(dm, 2, vertices, &joinSize, &join);CHKERRQ(ierr); 481 } 482 } 483 } 484 } 485 /* Create new pointSF from hashed claims */ 486 ierr = PetscHMapIGetSize(claimshash, &numLocalNew);CHKERRQ(ierr); 487 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 488 ierr = PetscMalloc1(numLeaves + numLocalNew, &localPointsNew);CHKERRQ(ierr); 489 ierr = PetscMalloc1(numLeaves + numLocalNew, &remotePointsNew);CHKERRQ(ierr); 490 for (p = 0; p < numLeaves; ++p) { 491 localPointsNew[p] = localPoints[p]; 492 remotePointsNew[p].index = remotePoints[p].index; 493 remotePointsNew[p].rank = remotePoints[p].rank; 494 } 495 p = numLeaves; 496 ierr = PetscHMapIGetKeys(claimshash, &p, localPointsNew);CHKERRQ(ierr); 497 ierr = PetscSortInt(numLocalNew, &localPointsNew[numLeaves]);CHKERRQ(ierr); 498 for (p = numLeaves; p < numLeaves + numLocalNew; ++p) { 499 ierr = PetscHMapIGet(claimshash, localPointsNew[p], &offset);CHKERRQ(ierr); 500 remotePointsNew[p] = claims[offset]; 501 } 502 ierr = PetscSFCreate(PetscObjectComm((PetscObject) dm), &sfPointNew);CHKERRQ(ierr); 503 ierr = PetscSFSetGraph(sfPointNew, pEnd-pStart, numLeaves+numLocalNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 504 ierr = DMSetPointSF(dm, sfPointNew);CHKERRQ(ierr); 505 ierr = PetscSFDestroy(&sfPointNew);CHKERRQ(ierr); 506 ierr = PetscHMapIDestroy(&claimshash);CHKERRQ(ierr); 507 } 508 ierr = PetscHMapIDestroy(&leafhash);CHKERRQ(ierr); 509 ierr = PetscHMapIJDestroy(&roothash);CHKERRQ(ierr); 510 ierr = PetscSectionDestroy(&candidateSection);CHKERRQ(ierr); 511 ierr = PetscSectionDestroy(&candidateSectionRemote);CHKERRQ(ierr); 512 ierr = PetscSectionDestroy(&claimSection);CHKERRQ(ierr); 513 ierr = PetscFree(candidates);CHKERRQ(ierr); 514 ierr = PetscFree(candidatesRemote);CHKERRQ(ierr); 515 ierr = PetscFree(claims);CHKERRQ(ierr); 516 ierr = PetscLogEventEnd(DMPLEX_InterpolateSF,dm,0,0,0);CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519 520 /*@C 521 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 522 523 Collective on DM 524 525 Input Parameters: 526 + dm - The DMPlex object with only cells and vertices 527 - dmInt - The interpolated DM 528 529 Output Parameter: 530 . dmInt - The complete DMPlex object 531 532 Level: intermediate 533 534 Notes: 535 It does not copy over the coordinates. 536 537 .keywords: mesh 538 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 539 @*/ 540 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 541 { 542 DM idm, odm = dm; 543 PetscSF sfPoint; 544 PetscInt depth, dim, d; 545 const char *name; 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 550 PetscValidPointer(dmInt, 2); 551 ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 552 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 553 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 554 if ((depth == dim) || (dim <= 1)) { 555 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 556 idm = dm; 557 } else { 558 for (d = 1; d < dim; ++d) { 559 /* Create interpolated mesh */ 560 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 561 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 562 ierr = DMSetDimension(idm, dim);CHKERRQ(ierr); 563 if (depth > 0) { 564 ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr); 565 ierr = DMGetPointSF(odm, &sfPoint);CHKERRQ(ierr); 566 ierr = DMPlexInterpolatePointSF(idm, sfPoint, depth);CHKERRQ(ierr); 567 } 568 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 569 odm = idm; 570 } 571 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 572 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 573 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 574 ierr = DMCopyLabels(dm, idm);CHKERRQ(ierr); 575 } 576 { 577 PetscBool isper; 578 const PetscReal *maxCell, *L; 579 const DMBoundaryType *bd; 580 581 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 582 ierr = DMSetPeriodicity(idm,isper,maxCell,L,bd);CHKERRQ(ierr); 583 } 584 *dmInt = idm; 585 ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr); 586 PetscFunctionReturn(0); 587 } 588 589 /*@ 590 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 591 592 Collective on DM 593 594 Input Parameter: 595 . dmA - The DMPlex object with initial coordinates 596 597 Output Parameter: 598 . dmB - The DMPlex object with copied coordinates 599 600 Level: intermediate 601 602 Note: This is typically used when adding pieces other than vertices to a mesh 603 604 .keywords: mesh 605 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection() 606 @*/ 607 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 608 { 609 Vec coordinatesA, coordinatesB; 610 VecType vtype; 611 PetscSection coordSectionA, coordSectionB; 612 PetscScalar *coordsA, *coordsB; 613 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 614 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE; 615 PetscBool lc = PETSC_FALSE; 616 PetscErrorCode ierr; 617 618 PetscFunctionBegin; 619 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 620 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 621 if (dmA == dmB) PetscFunctionReturn(0); 622 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 623 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 624 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); 625 ierr = DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA);CHKERRQ(ierr); 626 ierr = DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB);CHKERRQ(ierr); 627 ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 628 ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 629 if (coordSectionA == coordSectionB) PetscFunctionReturn(0); 630 ierr = PetscSectionGetNumFields(coordSectionA, &Nf);CHKERRQ(ierr); 631 if (!Nf) PetscFunctionReturn(0); 632 if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %D", Nf); 633 if (!coordSectionB) { 634 PetscInt dim; 635 636 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) coordSectionA), &coordSectionB);CHKERRQ(ierr); 637 ierr = DMGetCoordinateDim(dmA, &dim);CHKERRQ(ierr); 638 ierr = DMSetCoordinateSection(dmB, dim, coordSectionB);CHKERRQ(ierr); 639 ierr = PetscObjectDereference((PetscObject) coordSectionB);CHKERRQ(ierr); 640 } 641 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 642 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 643 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 644 ierr = PetscSectionGetChart(coordSectionA, &cS, &cE);CHKERRQ(ierr); 645 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 646 if ((cEndA-cStartA) != (cEndB-cStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cellls in first DM %d != %d in the second DM", cEndA-cStartA, cEndB-cStartB); 647 cS = cS - cStartA + cStartB; 648 cE = vEndB; 649 lc = PETSC_TRUE; 650 } else { 651 cS = vStartB; 652 cE = vEndB; 653 } 654 ierr = PetscSectionSetChart(coordSectionB, cS, cE);CHKERRQ(ierr); 655 for (v = vStartB; v < vEndB; ++v) { 656 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 657 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 658 } 659 if (lc) { /* localized coordinates */ 660 PetscInt c; 661 662 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 663 PetscInt dof; 664 665 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 666 ierr = PetscSectionSetDof(coordSectionB, c + cStartB, dof);CHKERRQ(ierr); 667 ierr = PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof);CHKERRQ(ierr); 668 } 669 } 670 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 671 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 672 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 673 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 674 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 675 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 676 ierr = VecGetBlockSize(coordinatesA, &d);CHKERRQ(ierr); 677 ierr = VecSetBlockSize(coordinatesB, d);CHKERRQ(ierr); 678 ierr = VecGetType(coordinatesA, &vtype);CHKERRQ(ierr); 679 ierr = VecSetType(coordinatesB, vtype);CHKERRQ(ierr); 680 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 681 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 682 for (v = 0; v < vEndB-vStartB; ++v) { 683 PetscInt offA, offB; 684 685 ierr = PetscSectionGetOffset(coordSectionA, v + vStartA, &offA);CHKERRQ(ierr); 686 ierr = PetscSectionGetOffset(coordSectionB, v + vStartB, &offB);CHKERRQ(ierr); 687 for (d = 0; d < spaceDim; ++d) { 688 coordsB[offB+d] = coordsA[offA+d]; 689 } 690 } 691 if (lc) { /* localized coordinates */ 692 PetscInt c; 693 694 for (c = cS-cStartB; c < cEndB-cStartB; c++) { 695 PetscInt dof, offA, offB; 696 697 ierr = PetscSectionGetOffset(coordSectionA, c + cStartA, &offA);CHKERRQ(ierr); 698 ierr = PetscSectionGetOffset(coordSectionB, c + cStartB, &offB);CHKERRQ(ierr); 699 ierr = PetscSectionGetDof(coordSectionA, c + cStartA, &dof);CHKERRQ(ierr); 700 ierr = PetscMemcpy(coordsB + offB,coordsA + offA,dof*sizeof(*coordsB));CHKERRQ(ierr); 701 } 702 } 703 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 704 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 705 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 706 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 710 /*@ 711 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 712 713 Collective on DM 714 715 Input Parameter: 716 . dm - The complete DMPlex object 717 718 Output Parameter: 719 . dmUnint - The DMPlex object with only cells and vertices 720 721 Level: intermediate 722 723 Notes: 724 It does not copy over the coordinates. 725 726 .keywords: mesh 727 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList(), DMPlexCopyCoordinates() 728 @*/ 729 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 730 { 731 DM udm; 732 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 733 PetscErrorCode ierr; 734 735 PetscFunctionBegin; 736 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 737 PetscValidPointer(dmUnint, 2); 738 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 739 if (dim <= 1) { 740 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 741 *dmUnint = dm; 742 PetscFunctionReturn(0); 743 } 744 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 745 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 746 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr); 747 ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr); 748 ierr = DMSetDimension(udm, dim);CHKERRQ(ierr); 749 ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr); 750 for (c = cStart; c < cEnd; ++c) { 751 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 752 753 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 754 for (cl = 0; cl < closureSize*2; cl += 2) { 755 const PetscInt p = closure[cl]; 756 757 if ((p >= vStart) && (p < vEnd)) ++coneSize; 758 } 759 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 760 ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr); 761 maxConeSize = PetscMax(maxConeSize, coneSize); 762 } 763 ierr = DMSetUp(udm);CHKERRQ(ierr); 764 ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr); 765 for (c = cStart; c < cEnd; ++c) { 766 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 767 768 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 769 for (cl = 0; cl < closureSize*2; cl += 2) { 770 const PetscInt p = closure[cl]; 771 772 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 773 } 774 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 775 ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr); 776 } 777 ierr = PetscFree(cone);CHKERRQ(ierr); 778 ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr); 779 ierr = DMPlexStratify(udm);CHKERRQ(ierr); 780 /* Reduce SF */ 781 { 782 PetscSF sfPoint, sfPointUn; 783 const PetscSFNode *remotePoints; 784 const PetscInt *localPoints; 785 PetscSFNode *remotePointsUn; 786 PetscInt *localPointsUn; 787 PetscInt vEnd, numRoots, numLeaves, l; 788 PetscInt numLeavesUn = 0, n = 0; 789 PetscErrorCode ierr; 790 791 /* Get original SF information */ 792 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 793 ierr = DMGetPointSF(udm, &sfPointUn);CHKERRQ(ierr); 794 ierr = DMPlexGetDepthStratum(dm, 0, NULL, &vEnd);CHKERRQ(ierr); 795 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 796 /* Allocate space for cells and vertices */ 797 for (l = 0; l < numLeaves; ++l) if (localPoints[l] < vEnd) numLeavesUn++; 798 /* Fill in leaves */ 799 if (vEnd >= 0) { 800 ierr = PetscMalloc1(numLeavesUn, &remotePointsUn);CHKERRQ(ierr); 801 ierr = PetscMalloc1(numLeavesUn, &localPointsUn);CHKERRQ(ierr); 802 for (l = 0; l < numLeaves; l++) { 803 if (localPoints[l] < vEnd) { 804 localPointsUn[n] = localPoints[l]; 805 remotePointsUn[n].rank = remotePoints[l].rank; 806 remotePointsUn[n].index = remotePoints[l].index; 807 ++n; 808 } 809 } 810 if (n != numLeavesUn) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %d != %d", n, numLeavesUn); 811 ierr = PetscSFSetGraph(sfPointUn, vEnd, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER);CHKERRQ(ierr); 812 } 813 } 814 { 815 PetscBool isper; 816 const PetscReal *maxCell, *L; 817 const DMBoundaryType *bd; 818 819 ierr = DMGetPeriodicity(dm,&isper,&maxCell,&L,&bd);CHKERRQ(ierr); 820 ierr = DMSetPeriodicity(udm,isper,maxCell,L,bd);CHKERRQ(ierr); 821 } 822 823 *dmUnint = udm; 824 PetscFunctionReturn(0); 825 } 826