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