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" 6 /* 7 DMPlexGetFaces - 8 9 Note: This will only work for cell-vertex meshes. 10 */ 11 static PetscErrorCode DMPlexGetFaces(DM dm, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 12 { 13 DM_Plex *mesh = (DM_Plex*) dm->data; 14 const PetscInt *cone = NULL; 15 PetscInt depth = 0, dim, coneSize; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 20 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 21 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 22 if (depth > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Faces can only be returned for cell-vertex meshes."); 23 if (!mesh->facesTmp) {ierr = PetscMalloc(PetscSqr(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)) * sizeof(PetscInt), &mesh->facesTmp);CHKERRQ(ierr);} 24 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 25 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 26 switch (dim) { 27 case 2: 28 switch (coneSize) { 29 case 3: 30 mesh->facesTmp[0] = cone[0]; mesh->facesTmp[1] = cone[1]; 31 mesh->facesTmp[2] = cone[1]; mesh->facesTmp[3] = cone[2]; 32 mesh->facesTmp[4] = cone[2]; mesh->facesTmp[5] = cone[0]; 33 *numFaces = 3; 34 *faceSize = 2; 35 *faces = mesh->facesTmp; 36 break; 37 case 4: 38 mesh->facesTmp[0] = cone[0]; mesh->facesTmp[1] = cone[1]; 39 mesh->facesTmp[2] = cone[1]; mesh->facesTmp[3] = cone[2]; 40 mesh->facesTmp[4] = cone[2]; mesh->facesTmp[5] = cone[3]; 41 mesh->facesTmp[6] = cone[3]; mesh->facesTmp[7] = cone[0]; 42 *numFaces = 4; 43 *faceSize = 2; 44 *faces = mesh->facesTmp; 45 break; 46 default: 47 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 48 } 49 break; 50 case 3: 51 switch (coneSize) { 52 case 3: 53 mesh->facesTmp[0] = cone[0]; mesh->facesTmp[1] = cone[1]; 54 mesh->facesTmp[2] = cone[1]; mesh->facesTmp[3] = cone[2]; 55 mesh->facesTmp[4] = cone[2]; mesh->facesTmp[5] = cone[0]; 56 *numFaces = 3; 57 *faceSize = 2; 58 *faces = mesh->facesTmp; 59 break; 60 case 4: 61 mesh->facesTmp[0] = cone[0]; mesh->facesTmp[1] = cone[1]; mesh->facesTmp[2] = cone[2]; 62 mesh->facesTmp[3] = cone[0]; mesh->facesTmp[4] = cone[2]; mesh->facesTmp[5] = cone[3]; 63 mesh->facesTmp[6] = cone[0]; mesh->facesTmp[7] = cone[3]; mesh->facesTmp[8] = cone[1]; 64 mesh->facesTmp[9] = cone[1]; mesh->facesTmp[10] = cone[3]; mesh->facesTmp[11] = cone[2]; 65 *numFaces = 4; 66 *faceSize = 3; 67 *faces = mesh->facesTmp; 68 break; 69 default: 70 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 71 } 72 break; 73 default: 74 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 75 } 76 PetscFunctionReturn(0); 77 } 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "DMPlexInterpolate_2D" 81 static PetscErrorCode DMPlexInterpolate_2D(DM dm, DM *dmInt) 82 { 83 DM idm; 84 DM_Plex *mesh; 85 PetscHashIJ edgeTable; 86 PetscInt *off; 87 PetscInt dim, numCells, cStart, cEnd, c, numVertices, vStart, vEnd; 88 PetscInt numEdges, firstEdge, edge, e; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 93 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 94 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 95 numCells = cEnd - cStart; 96 numVertices = vEnd - vStart; 97 firstEdge = numCells + numVertices; 98 numEdges = 0; 99 /* Count edges using algorithm from CreateNeighborCSR */ 100 ierr = DMPlexCreateNeighborCSR(dm, NULL, &off, NULL);CHKERRQ(ierr); 101 if (off) { 102 PetscInt numCorners = 0; 103 104 numEdges = off[numCells]/2; 105 #if 0 106 /* Account for boundary edges: \sum_c 3 - neighbors = 3*numCells - totalNeighbors */ 107 numEdges += 3*numCells - off[numCells]; 108 #else 109 /* Account for boundary edges: \sum_c #faces - #neighbors = \sum_c #cellVertices - #neighbors = totalCorners - totalNeighbors */ 110 for (c = cStart; c < cEnd; ++c) { 111 PetscInt coneSize; 112 113 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 114 numCorners += coneSize; 115 } 116 numEdges += numCorners - off[numCells]; 117 #endif 118 } 119 #if 0 120 /* Check Euler characteristic V - E + F = 1 */ 121 if (numVertices && (numVertices-numEdges+numCells != 1)) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Euler characteristic of mesh is %d != 1", numVertices-numEdges+numCells); 122 #endif 123 /* Create interpolated mesh */ 124 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 125 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 126 ierr = DMPlexSetDimension(idm, dim);CHKERRQ(ierr); 127 ierr = DMPlexSetChart(idm, 0, numCells+numVertices+numEdges);CHKERRQ(ierr); 128 for (c = 0; c < numCells; ++c) { 129 PetscInt numCorners; 130 131 ierr = DMPlexGetConeSize(dm, c, &numCorners);CHKERRQ(ierr); 132 ierr = DMPlexSetConeSize(idm, c, numCorners);CHKERRQ(ierr); 133 } 134 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 135 ierr = DMPlexSetConeSize(idm, e, 2);CHKERRQ(ierr); 136 } 137 ierr = DMSetUp(idm);CHKERRQ(ierr); 138 /* Get edge cones from subsets of cell vertices */ 139 ierr = PetscHashIJCreate(&edgeTable);CHKERRQ(ierr); 140 ierr = PetscHashIJSetMultivalued(edgeTable, PETSC_FALSE);CHKERRQ(ierr); 141 142 for (c = 0, edge = firstEdge; c < numCells; ++c) { 143 const PetscInt *cellFaces; 144 PetscInt numCellFaces, faceSize, cf; 145 146 ierr = DMPlexGetFaces(dm, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 147 if (faceSize != 2) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize); 148 for (cf = 0; cf < numCellFaces; ++cf) { 149 #if 1 150 PetscHashIJKey key; 151 152 key.i = PetscMin(cellFaces[cf*faceSize+0], cellFaces[cf*faceSize+1]); 153 key.j = PetscMax(cellFaces[cf*faceSize+0], cellFaces[cf*faceSize+1]); 154 ierr = PetscHashIJGet(edgeTable, key, &e);CHKERRQ(ierr); 155 if (e < 0) { 156 ierr = DMPlexSetCone(idm, edge, &cellFaces[cf*faceSize]);CHKERRQ(ierr); 157 ierr = PetscHashIJAdd(edgeTable, key, edge);CHKERRQ(ierr); 158 e = edge++; 159 } 160 #else 161 PetscBool found = PETSC_FALSE; 162 163 /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */ 164 for (e = firstEdge; e < edge; ++e) { 165 const PetscInt *cone; 166 167 ierr = DMPlexGetCone(idm, e, &cone);CHKERRQ(ierr); 168 if (((cellFaces[cf*faceSize+0] == cone[0]) && (cellFaces[cf*faceSize+1] == cone[1])) || 169 ((cellFaces[cf*faceSize+0] == cone[1]) && (cellFaces[cf*faceSize+1] == cone[0]))) { 170 found = PETSC_TRUE; 171 break; 172 } 173 } 174 if (!found) { 175 ierr = DMPlexSetCone(idm, edge, &cellFaces[cf*faceSize]);CHKERRQ(ierr); 176 ++edge; 177 } 178 #endif 179 ierr = DMPlexInsertCone(idm, c, cf, e);CHKERRQ(ierr); 180 } 181 } 182 if (edge != firstEdge+numEdges) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges); 183 ierr = PetscHashIJDestroy(&edgeTable);CHKERRQ(ierr); 184 ierr = PetscFree(off);CHKERRQ(ierr); 185 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 186 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 187 mesh = (DM_Plex*) (idm)->data; 188 /* Orient edges */ 189 for (c = 0; c < numCells; ++c) { 190 const PetscInt *cone = NULL, *cellFaces; 191 PetscInt coneSize, coff, numCellFaces, faceSize, cf; 192 193 ierr = DMPlexGetConeSize(idm, c, &coneSize);CHKERRQ(ierr); 194 ierr = DMPlexGetCone(idm, c, &cone);CHKERRQ(ierr); 195 ierr = PetscSectionGetOffset(mesh->coneSection, c, &coff);CHKERRQ(ierr); 196 ierr = DMPlexGetFaces(dm, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 197 if (coneSize != numCellFaces) SETERRQ3(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numCellFaces); 198 for (cf = 0; cf < numCellFaces; ++cf) { 199 const PetscInt *econe = NULL; 200 PetscInt esize; 201 202 ierr = DMPlexGetConeSize(idm, cone[cf], &esize);CHKERRQ(ierr); 203 ierr = DMPlexGetCone(idm, cone[cf], &econe);CHKERRQ(ierr); 204 if (esize != 2) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[cf]); 205 if ((cellFaces[cf*faceSize+0] == econe[0]) && (cellFaces[cf*faceSize+1] == econe[1])) { 206 /* Correctly oriented */ 207 mesh->coneOrientations[coff+cf] = 0; 208 } else if ((cellFaces[cf*faceSize+0] == econe[1]) && (cellFaces[cf*faceSize+1] == econe[0])) { 209 /* Start at index 1, and reverse orientation */ 210 mesh->coneOrientations[coff+cf] = -(1+1); 211 } 212 } 213 } 214 *dmInt = idm; 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "DMPlexInterpolate_3D" 220 static PetscErrorCode DMPlexInterpolate_3D(DM dm, DM *dmInt) 221 { 222 DM idm, fdm; 223 DM_Plex *mesh; 224 PetscInt *off; 225 const PetscInt numCorners = 4; 226 PetscInt dim, numCells, cStart, cEnd, c, numVertices, vStart, vEnd; 227 PetscInt numFaces, firstFace, face, f, numEdges, firstEdge, edge, e; 228 PetscErrorCode ierr; 229 230 PetscFunctionBegin; 231 { 232 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 233 ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 234 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 235 } 236 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 237 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 238 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 239 numCells = cEnd - cStart; 240 numVertices = vEnd - vStart; 241 firstFace = numCells + numVertices; 242 numFaces = 0; 243 /* Count faces using algorithm from CreateNeighborCSR */ 244 ierr = DMPlexCreateNeighborCSR(dm, NULL, &off, NULL);CHKERRQ(ierr); 245 if (off) { 246 numFaces = off[numCells]/2; 247 /* Account for boundary faces: \sum_c 4 - neighbors = 4*numCells - totalNeighbors */ 248 numFaces += 4*numCells - off[numCells]; 249 } 250 /* Use Euler characteristic to get edges V - E + F - C = 1 */ 251 firstEdge = firstFace + numFaces; 252 numEdges = numVertices + numFaces - numCells - 1; 253 /* Create interpolated mesh */ 254 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 255 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 256 ierr = DMPlexSetDimension(idm, dim);CHKERRQ(ierr); 257 ierr = DMPlexSetChart(idm, 0, numCells+numVertices+numFaces+numEdges);CHKERRQ(ierr); 258 for (c = 0; c < numCells; ++c) { 259 ierr = DMPlexSetConeSize(idm, c, numCorners);CHKERRQ(ierr); 260 } 261 for (f = firstFace; f < firstFace+numFaces; ++f) { 262 ierr = DMPlexSetConeSize(idm, f, 3);CHKERRQ(ierr); 263 } 264 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 265 ierr = DMPlexSetConeSize(idm, e, 2);CHKERRQ(ierr); 266 } 267 ierr = DMSetUp(idm);CHKERRQ(ierr); 268 /* Get face cones from subsets of cell vertices */ 269 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &fdm);CHKERRQ(ierr); 270 ierr = DMSetType(fdm, DMPLEX);CHKERRQ(ierr); 271 ierr = DMPlexSetDimension(fdm, dim);CHKERRQ(ierr); 272 ierr = DMPlexSetChart(fdm, numCells, firstFace+numFaces);CHKERRQ(ierr); 273 for (f = firstFace; f < firstFace+numFaces; ++f) { 274 ierr = DMPlexSetConeSize(fdm, f, 3);CHKERRQ(ierr); 275 } 276 ierr = DMSetUp(fdm);CHKERRQ(ierr); 277 for (c = 0, face = firstFace; c < numCells; ++c) { 278 const PetscInt *cellFaces; 279 PetscInt numCellFaces, faceSize, cf; 280 281 ierr = DMPlexGetFaces(dm, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 282 if (faceSize != 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Tetrahedra cannot have face of size %D", faceSize); 283 for (cf = 0; cf < numCellFaces; ++cf) { 284 PetscBool found = PETSC_FALSE; 285 286 /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */ 287 for (f = firstFace; f < face; ++f) { 288 const PetscInt *cone = NULL; 289 290 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 291 if (((cellFaces[cf*faceSize+0] == cone[0]) && (cellFaces[cf*faceSize+1] == cone[1]) && (cellFaces[cf*faceSize+2] == cone[2])) || 292 ((cellFaces[cf*faceSize+0] == cone[1]) && (cellFaces[cf*faceSize+1] == cone[2]) && (cellFaces[cf*faceSize+2] == cone[0])) || 293 ((cellFaces[cf*faceSize+0] == cone[2]) && (cellFaces[cf*faceSize+1] == cone[0]) && (cellFaces[cf*faceSize+2] == cone[1])) || 294 ((cellFaces[cf*faceSize+0] == cone[0]) && (cellFaces[cf*faceSize+1] == cone[2]) && (cellFaces[cf*faceSize+2] == cone[1])) || 295 ((cellFaces[cf*faceSize+0] == cone[2]) && (cellFaces[cf*faceSize+1] == cone[1]) && (cellFaces[cf*faceSize+2] == cone[0])) || 296 ((cellFaces[cf*faceSize+0] == cone[1]) && (cellFaces[cf*faceSize+1] == cone[0]) && (cellFaces[cf*faceSize+2] == cone[2]))) { 297 found = PETSC_TRUE; 298 break; 299 } 300 } 301 if (!found) { 302 ierr = DMPlexSetCone(idm, face, &cellFaces[cf*faceSize]);CHKERRQ(ierr); 303 /* Save the vertices for orientation calculation */ 304 ierr = DMPlexSetCone(fdm, face, &cellFaces[cf*faceSize]);CHKERRQ(ierr); 305 ++face; 306 } 307 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 308 } 309 } 310 if (face != firstFace+numFaces) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-firstFace, numFaces); 311 /* Get edge cones from subsets of face vertices */ 312 for (f = firstFace, edge = firstEdge; f < firstFace+numFaces; ++f) { 313 const PetscInt *cellFaces; 314 PetscInt numCellFaces, faceSize, cf; 315 316 ierr = DMPlexGetFaces(idm, f, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 317 if (faceSize != 2) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize); 318 for (cf = 0; cf < numCellFaces; ++cf) { 319 PetscBool found = PETSC_FALSE; 320 321 /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */ 322 for (e = firstEdge; e < edge; ++e) { 323 const PetscInt *cone = NULL; 324 325 ierr = DMPlexGetCone(idm, e, &cone);CHKERRQ(ierr); 326 if (((cellFaces[cf*faceSize+0] == cone[0]) && (cellFaces[cf*faceSize+1] == cone[1])) || 327 ((cellFaces[cf*faceSize+0] == cone[1]) && (cellFaces[cf*faceSize+1] == cone[0]))) { 328 found = PETSC_TRUE; 329 break; 330 } 331 } 332 if (!found) { 333 ierr = DMPlexSetCone(idm, edge, &cellFaces[cf*faceSize]);CHKERRQ(ierr); 334 ++edge; 335 } 336 ierr = DMPlexInsertCone(idm, f, cf, e);CHKERRQ(ierr); 337 } 338 } 339 if (edge != firstEdge+numEdges) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges); 340 ierr = PetscFree(off);CHKERRQ(ierr); 341 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 342 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 343 mesh = (DM_Plex*) (idm)->data; 344 /* Orient edges */ 345 for (f = firstFace; f < firstFace+numFaces; ++f) { 346 const PetscInt *cone, *cellFaces; 347 PetscInt coneSize, coff, numCellFaces, faceSize, cf; 348 349 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 350 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 351 ierr = PetscSectionGetOffset(mesh->coneSection, f, &coff);CHKERRQ(ierr); 352 ierr = DMPlexGetFaces(fdm, f, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 353 if (coneSize != numCellFaces) SETERRQ3(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid number of edges %D for face %D should be %D", coneSize, f, numCellFaces); 354 for (cf = 0; cf < numCellFaces; ++cf) { 355 const PetscInt *econe; 356 PetscInt esize; 357 358 ierr = DMPlexGetConeSize(idm, cone[cf], &esize);CHKERRQ(ierr); 359 ierr = DMPlexGetCone(idm, cone[cf], &econe);CHKERRQ(ierr); 360 if (esize != 2) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[cf]); 361 if ((cellFaces[cf*faceSize+0] == econe[0]) && (cellFaces[cf*faceSize+1] == econe[1])) { 362 /* Correctly oriented */ 363 mesh->coneOrientations[coff+cf] = 0; 364 } else if ((cellFaces[cf*faceSize+0] == econe[1]) && (cellFaces[cf*faceSize+1] == econe[0])) { 365 /* Start at index 1, and reverse orientation */ 366 mesh->coneOrientations[coff+cf] = -(1+1); 367 } 368 } 369 } 370 ierr = DMDestroy(&fdm);CHKERRQ(ierr); 371 /* Orient faces */ 372 for (c = 0; c < numCells; ++c) { 373 const PetscInt *cone, *cellFaces; 374 PetscInt coneSize, coff, numCellFaces, faceSize, cf; 375 376 ierr = DMPlexGetConeSize(idm, c, &coneSize);CHKERRQ(ierr); 377 ierr = DMPlexGetCone(idm, c, &cone);CHKERRQ(ierr); 378 ierr = PetscSectionGetOffset(mesh->coneSection, c, &coff);CHKERRQ(ierr); 379 ierr = DMPlexGetFaces(dm, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 380 if (coneSize != numCellFaces) SETERRQ3(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numCellFaces); 381 for (cf = 0; cf < numCellFaces; ++cf) { 382 PetscInt *origClosure = NULL, *closure; 383 PetscInt closureSize, i; 384 385 ierr = DMPlexGetTransitiveClosure(idm, cone[cf], PETSC_TRUE, &closureSize, &origClosure);CHKERRQ(ierr); 386 if (closureSize != 7) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid closure size %D for face %D should be 7", closureSize, cone[cf]); 387 for (i = 4; i < 7; ++i) { 388 if ((origClosure[i*2] < vStart) || (origClosure[i*2] >= vEnd)) SETERRQ3(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Invalid closure point %D should be a vertex in [%D, %D)", origClosure[i*2], vStart, vEnd); 389 } 390 closure = &origClosure[4*2]; 391 /* Remember that this is the orientation for edges, not vertices */ 392 if ((cellFaces[cf*faceSize+0] == closure[0*2]) && (cellFaces[cf*faceSize+1] == closure[1*2]) && (cellFaces[cf*faceSize+2] == closure[2*2])) { 393 /* Correctly oriented */ 394 mesh->coneOrientations[coff+cf] = 0; 395 } else if ((cellFaces[cf*faceSize+0] == closure[1*2]) && (cellFaces[cf*faceSize+1] == closure[2*2]) && (cellFaces[cf*faceSize+2] == closure[0*2])) { 396 /* Shifted by 1 */ 397 mesh->coneOrientations[coff+cf] = 1; 398 } else if ((cellFaces[cf*faceSize+0] == closure[2*2]) && (cellFaces[cf*faceSize+1] == closure[0*2]) && (cellFaces[cf*faceSize+2] == closure[1*2])) { 399 /* Shifted by 2 */ 400 mesh->coneOrientations[coff+cf] = 2; 401 } else if ((cellFaces[cf*faceSize+0] == closure[2*2]) && (cellFaces[cf*faceSize+1] == closure[1*2]) && (cellFaces[cf*faceSize+2] == closure[0*2])) { 402 /* Start at edge 1, and reverse orientation */ 403 mesh->coneOrientations[coff+cf] = -(1+1); 404 } else if ((cellFaces[cf*faceSize+0] == closure[1*2]) && (cellFaces[cf*faceSize+1] == closure[0*2]) && (cellFaces[cf*faceSize+2] == closure[2*2])) { 405 /* Start at index 0, and reverse orientation */ 406 mesh->coneOrientations[coff+cf] = -(0+1); 407 } else if ((cellFaces[cf*faceSize+0] == closure[0*2]) && (cellFaces[cf*faceSize+1] == closure[2*2]) && (cellFaces[cf*faceSize+2] == closure[1*2])) { 408 /* Start at index 2, and reverse orientation */ 409 mesh->coneOrientations[coff+cf] = -(2+1); 410 } else SETERRQ3(PetscObjectComm((PetscObject)idm), PETSC_ERR_PLIB, "Face %D did not match local face %D in cell %D for any orientation", cone[cf], cf, c); 411 ierr = DMPlexRestoreTransitiveClosure(idm, cone[cf], PETSC_TRUE, &closureSize, &origClosure);CHKERRQ(ierr); 412 } 413 } 414 { 415 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 416 ierr = DMView(idm, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 417 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 418 } 419 *dmInt = idm; 420 PetscFunctionReturn(0); 421 } 422 423 #undef __FUNCT__ 424 #define __FUNCT__ "DMPlexInterpolate" 425 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 426 { 427 PetscInt dim; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 432 switch (dim) { 433 case 2: 434 ierr = DMPlexInterpolate_2D(dm, dmInt);CHKERRQ(ierr);break; 435 case 3: 436 ierr = DMPlexInterpolate_3D(dm, dmInt);CHKERRQ(ierr);break; 437 default: 438 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No mesh interpolation support for dimension %D", dim); 439 } 440 PetscFunctionReturn(0); 441 } 442 443 #undef __FUNCT__ 444 #define __FUNCT__ "DMPlexCopyCoordinates" 445 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 446 { 447 Vec coordinatesA, coordinatesB; 448 PetscSection coordSectionA, coordSectionB; 449 PetscScalar *coordsA, *coordsB; 450 PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 451 PetscErrorCode ierr; 452 453 PetscFunctionBegin; 454 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 455 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 456 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); 457 ierr = DMPlexGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 458 ierr = DMPlexGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 459 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 460 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 461 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 462 ierr = PetscSectionSetChart(coordSectionB, vStartB, vEndB);CHKERRQ(ierr); 463 for (v = vStartB; v < vEndB; ++v) { 464 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 465 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 466 } 467 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 468 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 469 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 470 ierr = VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);CHKERRQ(ierr); 471 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 472 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 473 ierr = VecSetFromOptions(coordinatesB);CHKERRQ(ierr); 474 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 475 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 476 for (v = 0; v < vEndB-vStartB; ++v) { 477 for (d = 0; d < spaceDim; ++d) { 478 coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d]; 479 } 480 } 481 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 482 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 483 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 484 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 485 PetscFunctionReturn(0); 486 } 487