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