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