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 static PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[]) 10 { 11 const PetscInt *cone = NULL; 12 PetscInt *facesTmp; 13 PetscInt maxConeSize, maxSupportSize, coneSize; 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 18 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 19 ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);CHKERRQ(ierr); 20 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 21 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 22 switch (dim) { 23 case 2: 24 switch (coneSize) { 25 case 3: 26 if (faces) { 27 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 28 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 29 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 30 *faces = facesTmp; 31 } 32 if (numFaces) *numFaces = 3; 33 if (faceSize) *faceSize = 2; 34 break; 35 case 4: 36 /* Vertices follow right hand rule */ 37 if (faces) { 38 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 39 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 40 facesTmp[4] = cone[2]; facesTmp[5] = cone[3]; 41 facesTmp[6] = cone[3]; facesTmp[7] = cone[0]; 42 *faces = facesTmp; 43 } 44 if (numFaces) *numFaces = 4; 45 if (faceSize) *faceSize = 2; 46 if (faces) *faces = facesTmp; 47 break; 48 default: 49 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 50 } 51 break; 52 case 3: 53 switch (coneSize) { 54 case 3: 55 if (faces) { 56 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; 57 facesTmp[2] = cone[1]; facesTmp[3] = cone[2]; 58 facesTmp[4] = cone[2]; facesTmp[5] = cone[0]; 59 *faces = facesTmp; 60 } 61 if (numFaces) *numFaces = 3; 62 if (faceSize) *faceSize = 2; 63 if (faces) *faces = facesTmp; 64 break; 65 case 4: 66 /* Vertices of first face follow right hand rule and normal points away from last vertex */ 67 if (faces) { 68 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; 69 facesTmp[3] = cone[0]; facesTmp[4] = cone[3]; facesTmp[5] = cone[1]; 70 facesTmp[6] = cone[0]; facesTmp[7] = cone[2]; facesTmp[8] = cone[3]; 71 facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3]; 72 *faces = facesTmp; 73 } 74 if (numFaces) *numFaces = 4; 75 if (faceSize) *faceSize = 3; 76 if (faces) *faces = facesTmp; 77 break; 78 case 8: 79 if (faces) { 80 facesTmp[0] = cone[0]; facesTmp[1] = cone[1]; facesTmp[2] = cone[2]; facesTmp[3] = cone[3]; 81 facesTmp[4] = cone[4]; facesTmp[5] = cone[5]; facesTmp[6] = cone[6]; facesTmp[7] = cone[7]; 82 facesTmp[8] = cone[0]; facesTmp[9] = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; 83 facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; 84 facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; 85 facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; 86 *faces = facesTmp; 87 } 88 if (numFaces) *numFaces = 6; 89 if (faceSize) *faceSize = 4; 90 break; 91 default: 92 SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim); 93 } 94 break; 95 default: 96 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 97 } 98 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, &facesTmp);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "DMPlexInterpolateFaces_Internal" 104 /* This interpolates faces for cells at some stratum */ 105 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm) 106 { 107 DMLabel subpointMap; 108 PetscHashIJKL faceTable; 109 PetscInt *pStart, *pEnd; 110 PetscInt cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr); 115 /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */ 116 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 117 if (subpointMap) ++cellDim; 118 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 119 ++depth; 120 ++cellDepth; 121 cellDim -= depth - cellDepth; 122 ierr = PetscMalloc2(depth+1,PetscInt,&pStart,depth+1,PetscInt,&pEnd);CHKERRQ(ierr); 123 for (d = depth-1; d >= faceDepth; --d) { 124 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);CHKERRQ(ierr); 125 } 126 ierr = DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);CHKERRQ(ierr); 127 pEnd[faceDepth] = pStart[faceDepth]; 128 for (d = faceDepth-1; d >= 0; --d) { 129 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 130 } 131 if (pEnd[cellDepth] > pStart[cellDepth]) {ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr);} 132 if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 133 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 134 ierr = PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);CHKERRQ(ierr); 135 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 136 const PetscInt *cellFaces; 137 PetscInt numCellFaces, faceSize, cf, f; 138 139 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 140 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 141 for (cf = 0; cf < numCellFaces; ++cf) { 142 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 143 PetscHashIJKLKey key; 144 145 if (faceSize == 2) { 146 key.i = PetscMin(cellFace[0], cellFace[1]); 147 key.j = PetscMax(cellFace[0], cellFace[1]); 148 key.k = 0; 149 key.l = 0; 150 } else { 151 key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0; 152 ierr = PetscSortInt(faceSize, (PetscInt *) &key); 153 } 154 ierr = PetscHashIJKLGet(faceTable, key, &f);CHKERRQ(ierr); 155 if (f < 0) { 156 ierr = PetscHashIJKLAdd(faceTable, key, face);CHKERRQ(ierr); 157 f = face++; 158 } 159 } 160 } 161 pEnd[faceDepth] = face; 162 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 163 /* Count new points */ 164 for (d = 0; d <= depth; ++d) { 165 numPoints += pEnd[d]-pStart[d]; 166 } 167 ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr); 168 /* Set cone sizes */ 169 for (d = 0; d <= depth; ++d) { 170 PetscInt coneSize, p; 171 172 if (d == faceDepth) { 173 for (p = pStart[d]; p < pEnd[d]; ++p) { 174 /* I see no way to do this if we admit faces of different shapes */ 175 ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr); 176 } 177 } else if (d == cellDepth) { 178 for (p = pStart[d]; p < pEnd[d]; ++p) { 179 /* Number of cell faces may be different from number of cell vertices*/ 180 ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr); 181 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 182 } 183 } else { 184 for (p = pStart[d]; p < pEnd[d]; ++p) { 185 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 186 ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr); 187 } 188 } 189 } 190 ierr = DMSetUp(idm);CHKERRQ(ierr); 191 /* Get face cones from subsets of cell vertices */ 192 if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll); 193 ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr); 194 ierr = PetscHashIJKLSetMultivalued(faceTable, PETSC_FALSE);CHKERRQ(ierr); 195 for (d = depth; d > cellDepth; --d) { 196 const PetscInt *cone; 197 PetscInt p; 198 199 for (p = pStart[d]; p < pEnd[d]; ++p) { 200 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 201 ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr); 202 ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr); 203 ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr); 204 } 205 } 206 for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) { 207 const PetscInt *cellFaces; 208 PetscInt numCellFaces, faceSize, cf, f; 209 210 ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr); 211 if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll); 212 for (cf = 0; cf < numCellFaces; ++cf) { 213 const PetscInt *cellFace = &cellFaces[cf*faceSize]; 214 PetscHashIJKLKey key; 215 216 if (faceSize == 2) { 217 key.i = PetscMin(cellFace[0], cellFace[1]); 218 key.j = PetscMax(cellFace[0], cellFace[1]); 219 key.k = 0; 220 key.l = 0; 221 } else { 222 key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0; 223 ierr = PetscSortInt(faceSize, (PetscInt *) &key); 224 } 225 ierr = PetscHashIJKLGet(faceTable, key, &f);CHKERRQ(ierr); 226 if (f < 0) { 227 ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr); 228 ierr = PetscHashIJKLAdd(faceTable, key, face);CHKERRQ(ierr); 229 f = face++; 230 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 231 } else { 232 const PetscInt *cone; 233 PetscInt coneSize, ornt, i, j; 234 235 ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr); 236 /* Orient face: Do not allow reverse orientation at the first vertex */ 237 ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr); 238 ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr); 239 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); 240 /* - First find the initial vertex */ 241 for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break; 242 /* - Try forward comparison */ 243 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break; 244 if (j == faceSize) { 245 if ((faceSize == 2) && (i == 1)) ornt = -2; 246 else ornt = i; 247 } else { 248 /* - Try backward comparison */ 249 for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break; 250 if (j == faceSize) { 251 if (i == 0) ornt = -faceSize; 252 else ornt = -(i+1); 253 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation"); 254 } 255 ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr); 256 } 257 } 258 } 259 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]); 260 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 261 ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr); 262 ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr); 263 ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr); 264 ierr = DMPlexStratify(idm);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "DMPlexInterpolate" 270 /*@ 271 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 272 273 Collective on DM 274 275 Input Parameter: 276 . dmA - The DMPlex object with only cells and vertices 277 278 Output Parameter: 279 . dmB - The complete DMPlex object 280 281 Level: intermediate 282 283 .keywords: mesh 284 .seealso: DMPlexCreateFromCellList() 285 @*/ 286 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 287 { 288 DM idm, odm = dm; 289 PetscInt depth, dim, d; 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 294 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 295 if (dim <= 1) { 296 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 297 idm = dm; 298 } 299 for (d = 1; d < dim; ++d) { 300 /* Create interpolated mesh */ 301 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr); 302 ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr); 303 ierr = DMPlexSetDimension(idm, dim);CHKERRQ(ierr); 304 if (depth > 0) {ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr);} 305 if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 306 odm = idm; 307 } 308 *dmInt = idm; 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "DMPlexCopyCoordinates" 314 /*@ 315 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 316 317 Collective on DM 318 319 Input Parameter: 320 . dmA - The DMPlex object with initial coordinates 321 322 Output Parameter: 323 . dmB - The DMPlex object with copied coordinates 324 325 Level: intermediate 326 327 Note: This is typically used when adding pieces other than vertices to a mesh 328 329 .keywords: mesh 330 .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection() 331 @*/ 332 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 333 { 334 Vec coordinatesA, coordinatesB; 335 PetscSection coordSectionA, coordSectionB; 336 PetscScalar *coordsA, *coordsB; 337 PetscInt spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 if (dmA == dmB) PetscFunctionReturn(0); 342 ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr); 343 ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr); 344 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); 345 ierr = DMPlexGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr); 346 ierr = DMPlexGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr); 347 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 348 ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr); 349 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr); 350 ierr = PetscSectionSetChart(coordSectionB, vStartB, vEndB);CHKERRQ(ierr); 351 for (v = vStartB; v < vEndB; ++v) { 352 ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr); 353 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr); 354 } 355 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 356 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr); 357 ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr); 358 ierr = VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);CHKERRQ(ierr); 359 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 360 ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr); 361 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 362 ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr); 363 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 364 for (v = 0; v < vEndB-vStartB; ++v) { 365 for (d = 0; d < spaceDim; ++d) { 366 coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d]; 367 } 368 } 369 ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr); 370 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 371 ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr); 372 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375