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