1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /*I "petscdmfield.h" I*/ 3 #include <petscfe.h> 4 #include <petscdmplex.h> 5 #include <petscds.h> 6 7 typedef struct _n_DMField_DS { 8 PetscBool multifieldVec; 9 PetscInt height; /* Point height at which we want values and number of discretizations */ 10 PetscInt fieldNum; /* Number in DS of field which we evaluate */ 11 PetscObject *disc; /* Discretizations of this field at each height */ 12 Vec vec; /* Field values */ 13 DM dmDG; /* DM for the DG values */ 14 PetscObject *discDG; /* DG Discretizations of this field at each height */ 15 Vec vecDG; /* DG Field values */ 16 } DMField_DS; 17 18 static PetscErrorCode DMFieldDestroy_DS(DMField field) 19 { 20 DMField_DS *dsfield = (DMField_DS *)field->data; 21 PetscInt i; 22 23 PetscFunctionBegin; 24 PetscCall(VecDestroy(&dsfield->vec)); 25 for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->disc[i])); 26 PetscCall(PetscFree(dsfield->disc)); 27 PetscCall(VecDestroy(&dsfield->vecDG)); 28 if (dsfield->discDG) 29 for (i = 0; i < dsfield->height; i++) PetscCall(PetscObjectDereference(dsfield->discDG[i])); 30 PetscCall(PetscFree(dsfield->discDG)); 31 PetscCall(PetscFree(dsfield)); 32 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 static PetscErrorCode DMFieldView_DS(DMField field, PetscViewer viewer) 36 { 37 DMField_DS *dsfield = (DMField_DS *)field->data; 38 PetscObject disc; 39 PetscBool iascii; 40 41 PetscFunctionBegin; 42 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 43 disc = dsfield->disc[0]; 44 if (iascii) { 45 PetscCall(PetscViewerASCIIPrintf(viewer, "PetscDS field %" PetscInt_FMT "\n", dsfield->fieldNum)); 46 PetscCall(PetscViewerASCIIPushTab(viewer)); 47 PetscCall(PetscObjectView(disc, viewer)); 48 PetscCall(PetscViewerASCIIPopTab(viewer)); 49 } 50 PetscCall(PetscViewerASCIIPushTab(viewer)); 51 PetscCheck(!dsfield->multifieldVec, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "View of subfield not implemented yet"); 52 PetscCall(VecView(dsfield->vec, viewer)); 53 PetscCall(PetscViewerASCIIPopTab(viewer)); 54 PetscFunctionReturn(PETSC_SUCCESS); 55 } 56 57 static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject discList[], PetscObject *disc) 58 { 59 PetscFunctionBegin; 60 if (!discList[height]) { 61 PetscClassId id; 62 63 PetscCall(PetscObjectGetClassId(discList[0], &id)); 64 if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height])); 65 } 66 *disc = discList[height]; 67 PetscFunctionReturn(PETSC_SUCCESS); 68 } 69 70 /* y[m,c] = A[m,n,c] . b[n] */ 71 #define DMFieldDSdot(y, A, b, m, n, c, cast) \ 72 do { \ 73 PetscInt _i, _j, _k; \ 74 for (_i = 0; _i < (m); _i++) { \ 75 for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] = 0.; \ 76 for (_j = 0; _j < (n); _j++) { \ 77 for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \ 78 } \ 79 } \ 80 } while (0) 81 82 /* 83 Since this is used for coordinates, we need to allow for the possibility that values come from multiple sections/Vecs, so that we can have DG version of the coordinates for periodicity. This reproduces DMPlexGetCellCoordinates_Internal(). 84 */ 85 PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[]) 86 { 87 DMField_DS *dsfield = (DMField_DS *)field->data; 88 DM fdm = dsfield->dmDG; 89 PetscSection s = NULL; 90 const PetscScalar *cvalues; 91 PetscInt pStart, pEnd; 92 93 PetscFunctionBeginHot; 94 *isDG = PETSC_FALSE; 95 *Nc = 0; 96 *array = NULL; 97 *values = NULL; 98 /* Check for cellwise section */ 99 if (fdm) PetscCall(DMGetLocalSection(fdm, &s)); 100 if (!s) goto cg; 101 /* Check that the cell exists in the cellwise section */ 102 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 103 if (cell < pStart || cell >= pEnd) goto cg; 104 /* Check for cellwise coordinates for this cell */ 105 PetscCall(PetscSectionGetDof(s, cell, Nc)); 106 if (!*Nc) goto cg; 107 /* Check for cellwise coordinates */ 108 if (!dsfield->vecDG) goto cg; 109 /* Get cellwise coordinates */ 110 PetscCall(VecGetArrayRead(dsfield->vecDG, array)); 111 PetscCall(DMPlexPointLocalRead(fdm, cell, *array, &cvalues)); 112 PetscCall(DMGetWorkArray(fdm, *Nc, MPIU_SCALAR, values)); 113 PetscCall(PetscArraycpy(*values, cvalues, *Nc)); 114 PetscCall(VecRestoreArrayRead(dsfield->vecDG, array)); 115 *isDG = PETSC_TRUE; 116 PetscFunctionReturn(PETSC_SUCCESS); 117 cg: 118 /* Use continuous values */ 119 PetscCall(DMFieldGetDM(field, &fdm)); 120 PetscCall(DMGetLocalSection(fdm, &s)); 121 PetscCall(PetscSectionGetField(s, dsfield->fieldNum, &s)); 122 PetscCall(DMPlexVecGetClosure(fdm, s, dsfield->vec, cell, Nc, values)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 PetscErrorCode DMFieldRestoreClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[]) 127 { 128 DMField_DS *dsfield = (DMField_DS *)field->data; 129 DM fdm; 130 PetscSection s; 131 132 PetscFunctionBeginHot; 133 if (*isDG) { 134 PetscCall(DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values)); 135 } else { 136 PetscCall(DMFieldGetDM(field, &fdm)); 137 PetscCall(DMGetLocalSection(fdm, &s)); 138 PetscCall(DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, (PetscScalar **)values)); 139 } 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 143 /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */ 144 static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H) 145 { 146 DMField_DS *dsfield = (DMField_DS *)field->data; 147 DM dm; 148 PetscObject disc; 149 PetscClassId classid; 150 PetscInt nq, nc, dim, meshDim, numCells; 151 PetscSection section; 152 const PetscReal *qpoints; 153 PetscBool isStride; 154 const PetscInt *points = NULL; 155 PetscInt sfirst = -1, stride = -1; 156 157 PetscFunctionBeginHot; 158 dm = field->dm; 159 nc = field->numComponents; 160 PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL)); 161 PetscCall(DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc)); 162 PetscCall(DMGetDimension(dm, &meshDim)); 163 PetscCall(DMGetLocalSection(dm, §ion)); 164 PetscCall(PetscSectionGetField(section, dsfield->fieldNum, §ion)); 165 PetscCall(PetscObjectGetClassId(disc, &classid)); 166 /* TODO: batch */ 167 PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride)); 168 PetscCall(ISGetLocalSize(pointIS, &numCells)); 169 if (isStride) PetscCall(ISStrideGetInfo(pointIS, &sfirst, &stride)); 170 else PetscCall(ISGetIndices(pointIS, &points)); 171 if (classid == PETSCFE_CLASSID) { 172 PetscFE fe = (PetscFE)disc; 173 PetscInt feDim, i; 174 PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1)); 175 PetscTabulation T; 176 177 PetscCall(PetscFEGetDimension(fe, &feDim)); 178 PetscCall(PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T)); 179 for (i = 0; i < numCells; i++) { 180 PetscInt c = isStride ? (sfirst + i * stride) : points[i]; 181 PetscInt closureSize; 182 const PetscScalar *array; 183 PetscScalar *elem = NULL; 184 PetscBool isDG; 185 186 PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem)); 187 if (B) { 188 /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */ 189 if (type == PETSC_SCALAR) { 190 PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i]; 191 192 DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar)); 193 } else { 194 PetscReal *cB = &((PetscReal *)B)[nc * nq * i]; 195 196 DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart); 197 } 198 } 199 if (D) { 200 if (type == PETSC_SCALAR) { 201 PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i]; 202 203 DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar)); 204 } else { 205 PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i]; 206 207 DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart); 208 } 209 } 210 if (H) { 211 if (type == PETSC_SCALAR) { 212 PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i]; 213 214 DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar)); 215 } else { 216 PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i]; 217 218 DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart); 219 } 220 } 221 PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem)); 222 } 223 PetscCall(PetscTabulationDestroy(&T)); 224 } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented"); 225 if (!isStride) PetscCall(ISRestoreIndices(pointIS, &points)); 226 PetscFunctionReturn(PETSC_SUCCESS); 227 } 228 229 static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 230 { 231 DMField_DS *dsfield = (DMField_DS *)field->data; 232 PetscSF cellSF = NULL; 233 const PetscSFNode *cells; 234 PetscInt c, nFound, numCells, feDim, nc; 235 const PetscInt *cellDegrees; 236 const PetscScalar *pointsArray; 237 PetscScalar *cellPoints; 238 PetscInt gatherSize, gatherMax; 239 PetscInt dim, dimR, offset; 240 MPI_Datatype pointType; 241 PetscObject cellDisc; 242 PetscFE cellFE; 243 PetscClassId discID; 244 PetscReal *coordsReal, *coordsRef; 245 PetscSection section; 246 PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL; 247 PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL; 248 PetscReal *v, *J, *invJ, *detJ; 249 250 PetscFunctionBegin; 251 nc = field->numComponents; 252 PetscCall(DMGetLocalSection(field->dm, §ion)); 253 PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc)); 254 PetscCall(PetscObjectGetClassId(cellDisc, &discID)); 255 PetscCheck(discID == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization type not supported"); 256 cellFE = (PetscFE)cellDisc; 257 PetscCall(PetscFEGetDimension(cellFE, &feDim)); 258 PetscCall(DMGetCoordinateDim(field->dm, &dim)); 259 PetscCall(DMGetDimension(field->dm, &dimR)); 260 PetscCall(DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF)); 261 PetscCall(PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells)); 262 for (c = 0; c < nFound; c++) PetscCheck(cells[c].index >= 0, PetscObjectComm((PetscObject)points), PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " could not be located", c); 263 PetscCall(PetscSFComputeDegreeBegin(cellSF, &cellDegrees)); 264 PetscCall(PetscSFComputeDegreeEnd(cellSF, &cellDegrees)); 265 for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) { 266 gatherMax = PetscMax(gatherMax, cellDegrees[c]); 267 gatherSize += cellDegrees[c]; 268 } 269 PetscCall(PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef)); 270 PetscCall(PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ)); 271 if (datatype == PETSC_SCALAR) PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs)); 272 else PetscCall(PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr)); 273 274 PetscCallMPI(MPI_Type_contiguous(dim, MPIU_SCALAR, &pointType)); 275 PetscCallMPI(MPI_Type_commit(&pointType)); 276 PetscCall(VecGetArrayRead(points, &pointsArray)); 277 PetscCall(PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints)); 278 PetscCall(PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints)); 279 PetscCall(VecRestoreArrayRead(points, &pointsArray)); 280 for (c = 0, offset = 0; c < numCells; c++) { 281 PetscInt nq = cellDegrees[c], p; 282 283 if (nq) { 284 PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1)); 285 PetscTabulation T; 286 PetscQuadrature quad; 287 const PetscScalar *array; 288 PetscScalar *elem = NULL; 289 PetscReal *quadPoints; 290 PetscBool isDG; 291 PetscInt closureSize, d, e, f, g; 292 293 for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]); 294 PetscCall(DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef)); 295 PetscCall(PetscFECreateTabulation(cellFE, 1, nq, coordsRef, K, &T)); 296 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad)); 297 PetscCall(PetscMalloc1(dimR * nq, &quadPoints)); 298 for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p]; 299 PetscCall(PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL)); 300 PetscCall(DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ)); 301 PetscCall(PetscQuadratureDestroy(&quad)); 302 PetscCall(DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem)); 303 if (B) { 304 if (datatype == PETSC_SCALAR) { 305 PetscScalar *cB = &cellBs[nc * offset]; 306 307 DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar)); 308 } else { 309 PetscReal *cB = &cellBr[nc * offset]; 310 311 DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart); 312 } 313 } 314 if (D) { 315 if (datatype == PETSC_SCALAR) { 316 PetscScalar *cD = &cellDs[nc * dim * offset]; 317 318 DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar)); 319 for (p = 0; p < nq; p++) { 320 for (g = 0; g < nc; g++) { 321 PetscScalar vs[3]; 322 323 for (d = 0; d < dimR; d++) { 324 vs[d] = 0.; 325 for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e]; 326 } 327 for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = vs[d]; 328 } 329 } 330 } else { 331 PetscReal *cD = &cellDr[nc * dim * offset]; 332 333 DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart); 334 for (p = 0; p < nq; p++) { 335 for (g = 0; g < nc; g++) { 336 for (d = 0; d < dimR; d++) { 337 v[d] = 0.; 338 for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e]; 339 } 340 for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = v[d]; 341 } 342 } 343 } 344 } 345 if (H) { 346 if (datatype == PETSC_SCALAR) { 347 PetscScalar *cH = &cellHs[nc * dim * dim * offset]; 348 349 DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar)); 350 for (p = 0; p < nq; p++) { 351 for (g = 0; g < nc * dimR; g++) { 352 PetscScalar vs[3]; 353 354 for (d = 0; d < dimR; d++) { 355 vs[d] = 0.; 356 for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e]; 357 } 358 for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = vs[d]; 359 } 360 for (g = 0; g < nc; g++) { 361 for (f = 0; f < dimR; f++) { 362 PetscScalar vs[3]; 363 364 for (d = 0; d < dimR; d++) { 365 vs[d] = 0.; 366 for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f]; 367 } 368 for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d]; 369 } 370 } 371 } 372 } else { 373 PetscReal *cH = &cellHr[nc * dim * dim * offset]; 374 375 DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart); 376 for (p = 0; p < nq; p++) { 377 for (g = 0; g < nc * dimR; g++) { 378 for (d = 0; d < dimR; d++) { 379 v[d] = 0.; 380 for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e]; 381 } 382 for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = v[d]; 383 } 384 for (g = 0; g < nc; g++) { 385 for (f = 0; f < dimR; f++) { 386 for (d = 0; d < dimR; d++) { 387 v[d] = 0.; 388 for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f]; 389 } 390 for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = v[d]; 391 } 392 } 393 } 394 } 395 } 396 PetscCall(DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem)); 397 PetscCall(PetscTabulationDestroy(&T)); 398 } 399 offset += nq; 400 } 401 { 402 MPI_Datatype origtype; 403 404 if (datatype == PETSC_SCALAR) { 405 origtype = MPIU_SCALAR; 406 } else { 407 origtype = MPIU_REAL; 408 } 409 if (B) { 410 MPI_Datatype Btype; 411 412 PetscCallMPI(MPI_Type_contiguous(nc, origtype, &Btype)); 413 PetscCallMPI(MPI_Type_commit(&Btype)); 414 PetscCall(PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B)); 415 PetscCall(PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B)); 416 PetscCallMPI(MPI_Type_free(&Btype)); 417 } 418 if (D) { 419 MPI_Datatype Dtype; 420 421 PetscCallMPI(MPI_Type_contiguous(nc * dim, origtype, &Dtype)); 422 PetscCallMPI(MPI_Type_commit(&Dtype)); 423 PetscCall(PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D)); 424 PetscCall(PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D)); 425 PetscCallMPI(MPI_Type_free(&Dtype)); 426 } 427 if (H) { 428 MPI_Datatype Htype; 429 430 PetscCallMPI(MPI_Type_contiguous(nc * dim * dim, origtype, &Htype)); 431 PetscCallMPI(MPI_Type_commit(&Htype)); 432 PetscCall(PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H)); 433 PetscCall(PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H)); 434 PetscCallMPI(MPI_Type_free(&Htype)); 435 } 436 } 437 PetscCall(PetscFree4(v, J, invJ, detJ)); 438 PetscCall(PetscFree3(cellBr, cellDr, cellHr)); 439 PetscCall(PetscFree3(cellBs, cellDs, cellHs)); 440 PetscCall(PetscFree3(cellPoints, coordsReal, coordsRef)); 441 PetscCallMPI(MPI_Type_free(&pointType)); 442 PetscCall(PetscSFDestroy(&cellSF)); 443 PetscFunctionReturn(PETSC_SUCCESS); 444 } 445 446 static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H) 447 { 448 DMField_DS *dsfield = (DMField_DS *)field->data; 449 PetscInt h, imin; 450 PetscInt dim; 451 PetscClassId id; 452 PetscQuadrature quad = NULL; 453 PetscInt maxDegree; 454 PetscFEGeom *geom; 455 PetscInt Nq, Nc, dimC, qNc, N; 456 PetscInt numPoints; 457 void *qB = NULL, *qD = NULL, *qH = NULL; 458 const PetscReal *weights; 459 MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL; 460 PetscObject disc; 461 DMField coordField; 462 463 PetscFunctionBegin; 464 Nc = field->numComponents; 465 PetscCall(DMGetCoordinateDim(field->dm, &dimC)); 466 PetscCall(DMGetDimension(field->dm, &dim)); 467 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 468 PetscCall(ISGetMinMax(pointIS, &imin, NULL)); 469 for (h = 0; h < dsfield->height; h++) { 470 PetscInt hEnd; 471 472 PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd)); 473 if (imin < hEnd) break; 474 } 475 dim -= h; 476 PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc)); 477 PetscCall(PetscObjectGetClassId(disc, &id)); 478 PetscCheck(id == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported"); 479 PetscCall(DMGetCoordinateField(field->dm, &coordField)); 480 PetscCall(DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree)); 481 if (maxDegree <= 1) PetscCall(DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad)); 482 if (!quad) PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, &quad)); 483 PetscCheck(quad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages"); 484 PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom)); 485 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights)); 486 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components"); 487 N = numPoints * Nq * Nc; 488 if (B) PetscCall(DMGetWorkArray(field->dm, N, mpitype, &qB)); 489 if (D) PetscCall(DMGetWorkArray(field->dm, N * dimC, mpitype, &qD)); 490 if (H) PetscCall(DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH)); 491 PetscCall(DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH)); 492 if (B) { 493 PetscInt i, j, k; 494 495 if (type == PETSC_SCALAR) { 496 PetscScalar *sB = (PetscScalar *)B; 497 PetscScalar *sqB = (PetscScalar *)qB; 498 499 for (i = 0; i < numPoints; i++) { 500 PetscReal vol = 0.; 501 502 for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.; 503 for (k = 0; k < Nq; k++) { 504 vol += geom->detJ[i * Nq + k] * weights[k]; 505 for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j]; 506 } 507 for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol; 508 } 509 } else { 510 PetscReal *rB = (PetscReal *)B; 511 PetscReal *rqB = (PetscReal *)qB; 512 513 for (i = 0; i < numPoints; i++) { 514 PetscReal vol = 0.; 515 516 for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.; 517 for (k = 0; k < Nq; k++) { 518 vol += geom->detJ[i * Nq + k] * weights[k]; 519 for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j]; 520 } 521 for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol; 522 } 523 } 524 } 525 if (D) { 526 PetscInt i, j, k, l, m; 527 528 if (type == PETSC_SCALAR) { 529 PetscScalar *sD = (PetscScalar *)D; 530 PetscScalar *sqD = (PetscScalar *)qD; 531 532 for (i = 0; i < numPoints; i++) { 533 PetscReal vol = 0.; 534 535 for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.; 536 for (k = 0; k < Nq; k++) { 537 vol += geom->detJ[i * Nq + k] * weights[k]; 538 for (j = 0; j < Nc; j++) { 539 PetscScalar pD[3] = {0., 0., 0.}; 540 541 for (l = 0; l < dimC; l++) { 542 for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m]; 543 } 544 for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l]; 545 } 546 } 547 for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol; 548 } 549 } else { 550 PetscReal *rD = (PetscReal *)D; 551 PetscReal *rqD = (PetscReal *)qD; 552 553 for (i = 0; i < numPoints; i++) { 554 PetscReal vol = 0.; 555 556 for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.; 557 for (k = 0; k < Nq; k++) { 558 vol += geom->detJ[i * Nq + k] * weights[k]; 559 for (j = 0; j < Nc; j++) { 560 PetscReal pD[3] = {0., 0., 0.}; 561 562 for (l = 0; l < dimC; l++) { 563 for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m]; 564 } 565 for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l]; 566 } 567 } 568 for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol; 569 } 570 } 571 } 572 if (H) { 573 PetscInt i, j, k, l, m, q, r; 574 575 if (type == PETSC_SCALAR) { 576 PetscScalar *sH = (PetscScalar *)H; 577 PetscScalar *sqH = (PetscScalar *)qH; 578 579 for (i = 0; i < numPoints; i++) { 580 PetscReal vol = 0.; 581 582 for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.; 583 for (k = 0; k < Nq; k++) { 584 const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC]; 585 586 vol += geom->detJ[i * Nq + k] * weights[k]; 587 for (j = 0; j < Nc; j++) { 588 PetscScalar pH[3][3] = { 589 {0., 0., 0.}, 590 {0., 0., 0.}, 591 {0., 0., 0.} 592 }; 593 const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC]; 594 595 for (l = 0; l < dimC; l++) { 596 for (m = 0; m < dimC; m++) { 597 for (q = 0; q < dim; q++) { 598 for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r]; 599 } 600 } 601 } 602 for (l = 0; l < dimC; l++) { 603 for (m = 0; m < dimC; m++) sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m]; 604 } 605 } 606 } 607 for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol; 608 } 609 } else { 610 PetscReal *rH = (PetscReal *)H; 611 PetscReal *rqH = (PetscReal *)qH; 612 613 for (i = 0; i < numPoints; i++) { 614 PetscReal vol = 0.; 615 616 for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.; 617 for (k = 0; k < Nq; k++) { 618 const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC]; 619 620 vol += geom->detJ[i * Nq + k] * weights[k]; 621 for (j = 0; j < Nc; j++) { 622 PetscReal pH[3][3] = { 623 {0., 0., 0.}, 624 {0., 0., 0.}, 625 {0., 0., 0.} 626 }; 627 const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC]; 628 629 for (l = 0; l < dimC; l++) { 630 for (m = 0; m < dimC; m++) { 631 for (q = 0; q < dim; q++) { 632 for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r]; 633 } 634 } 635 } 636 for (l = 0; l < dimC; l++) { 637 for (m = 0; m < dimC; m++) rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m]; 638 } 639 } 640 } 641 for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol; 642 } 643 } 644 } 645 if (B) PetscCall(DMRestoreWorkArray(field->dm, N, mpitype, &qB)); 646 if (D) PetscCall(DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD)); 647 if (H) PetscCall(DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH)); 648 PetscCall(PetscFEGeomDestroy(&geom)); 649 PetscCall(PetscQuadratureDestroy(&quad)); 650 PetscFunctionReturn(PETSC_SUCCESS); 651 } 652 653 static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree) 654 { 655 DMField_DS *dsfield; 656 PetscObject disc; 657 PetscInt h, imin, imax; 658 PetscClassId id; 659 660 PetscFunctionBegin; 661 dsfield = (DMField_DS *)field->data; 662 PetscCall(ISGetMinMax(pointIS, &imin, &imax)); 663 if (imin >= imax) { 664 h = 0; 665 } else { 666 for (h = 0; h < dsfield->height; h++) { 667 PetscInt hEnd; 668 669 PetscCall(DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd)); 670 if (imin < hEnd) break; 671 } 672 } 673 PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc)); 674 PetscCall(PetscObjectGetClassId(disc, &id)); 675 if (id == PETSCFE_CLASSID) { 676 PetscFE fe = (PetscFE)disc; 677 PetscSpace sp; 678 679 PetscCall(PetscFEGetBasisSpace(fe, &sp)); 680 PetscCall(PetscSpaceGetDegree(sp, minDegree, maxDegree)); 681 } 682 PetscFunctionReturn(PETSC_SUCCESS); 683 } 684 685 PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad) 686 { 687 DM dm = field->dm; 688 const PetscInt *points; 689 DMPolytopeType ct; 690 PetscInt dim, n; 691 PetscBool isplex; 692 693 PetscFunctionBegin; 694 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex)); 695 PetscCall(ISGetLocalSize(pointIS, &n)); 696 if (isplex && n) { 697 PetscCall(DMGetDimension(dm, &dim)); 698 PetscCall(ISGetIndices(pointIS, &points)); 699 PetscCall(DMPlexGetCellType(dm, points[0], &ct)); 700 switch (ct) { 701 case DM_POLYTOPE_TRIANGLE: 702 case DM_POLYTOPE_TETRAHEDRON: 703 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad)); 704 break; 705 default: 706 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad)); 707 } 708 PetscCall(ISRestoreIndices(pointIS, &points)); 709 } else PetscCall(DMFieldCreateDefaultQuadrature(field, pointIS, quad)); 710 PetscFunctionReturn(PETSC_SUCCESS); 711 } 712 713 static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad) 714 { 715 PetscInt h, dim, imax, imin, cellHeight; 716 DM dm; 717 DMField_DS *dsfield; 718 PetscObject disc; 719 PetscFE fe; 720 PetscClassId id; 721 722 PetscFunctionBegin; 723 dm = field->dm; 724 dsfield = (DMField_DS *)field->data; 725 PetscCall(ISGetMinMax(pointIS, &imin, &imax)); 726 PetscCall(DMGetDimension(dm, &dim)); 727 for (h = 0; h <= dim; h++) { 728 PetscInt hStart, hEnd; 729 730 PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd)); 731 if (imax >= hStart && imin < hEnd) break; 732 } 733 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 734 h -= cellHeight; 735 *quad = NULL; 736 if (h < dsfield->height) { 737 PetscCall(DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc)); 738 PetscCall(PetscObjectGetClassId(disc, &id)); 739 if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS); 740 fe = (PetscFE)disc; 741 PetscCall(PetscFEGetQuadrature(fe, quad)); 742 PetscCall(PetscObjectReference((PetscObject)*quad)); 743 } 744 PetscFunctionReturn(PETSC_SUCCESS); 745 } 746 747 static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom) 748 { 749 const PetscInt *points; 750 PetscInt p, dim, dE, numFaces, Nq; 751 PetscInt maxDegree; 752 DMLabel depthLabel; 753 IS cellIS; 754 DM dm = field->dm; 755 756 PetscFunctionBegin; 757 dim = geom->dim; 758 dE = geom->dimEmbed; 759 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 760 PetscCall(DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS)); 761 PetscCall(DMFieldGetDegree(field, cellIS, NULL, &maxDegree)); 762 PetscCall(ISGetIndices(pointIS, &points)); 763 numFaces = geom->numCells; 764 Nq = geom->numPoints; 765 /* First, set local faces and flip normals so that they are outward for the first supporting cell */ 766 for (p = 0; p < numFaces; p++) { 767 PetscInt point = points[p]; 768 PetscInt suppSize, s, coneSize, c, numChildren; 769 const PetscInt *supp, *cone, *ornt; 770 771 PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL)); 772 PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 773 PetscCall(DMPlexGetSupportSize(dm, point, &suppSize)); 774 PetscCheck(suppSize <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, suppSize); 775 if (!suppSize) continue; 776 PetscCall(DMPlexGetSupport(dm, point, &supp)); 777 for (s = 0; s < suppSize; ++s) { 778 PetscCall(DMPlexGetConeSize(dm, supp[s], &coneSize)); 779 PetscCall(DMPlexGetCone(dm, supp[s], &cone)); 780 for (c = 0; c < coneSize; ++c) 781 if (cone[c] == point) break; 782 PetscCheck(c != coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid connectivity: point %" PetscInt_FMT " not found in cone of support point %" PetscInt_FMT, point, supp[s]); 783 geom->face[p][s] = c; 784 } 785 PetscCall(DMPlexGetConeOrientation(dm, supp[0], &ornt)); 786 if (ornt[geom->face[p][0]] < 0) { 787 PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d; 788 789 for (q = 0; q < Np; ++q) 790 for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d]; 791 } 792 } 793 if (maxDegree <= 1) { 794 PetscInt numCells, offset, *cells; 795 PetscFEGeom *cellGeom; 796 IS suppIS; 797 798 for (p = 0, numCells = 0; p < numFaces; p++) { 799 PetscInt point = points[p]; 800 PetscInt numSupp, numChildren; 801 802 PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL)); 803 PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 804 PetscCall(DMPlexGetSupportSize(dm, point, &numSupp)); 805 PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp); 806 numCells += numSupp; 807 } 808 PetscCall(PetscMalloc1(numCells, &cells)); 809 for (p = 0, offset = 0; p < numFaces; p++) { 810 PetscInt point = points[p]; 811 PetscInt numSupp, s; 812 const PetscInt *supp; 813 814 PetscCall(DMPlexGetSupportSize(dm, point, &numSupp)); 815 PetscCall(DMPlexGetSupport(dm, point, &supp)); 816 for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s]; 817 } 818 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS)); 819 PetscCall(DMFieldCreateFEGeom(field, suppIS, quad, PETSC_FALSE, &cellGeom)); 820 for (p = 0, offset = 0; p < numFaces; p++) { 821 PetscInt point = points[p]; 822 PetscInt numSupp, s, q; 823 const PetscInt *supp; 824 825 PetscCall(DMPlexGetSupportSize(dm, point, &numSupp)); 826 PetscCall(DMPlexGetSupport(dm, point, &supp)); 827 for (s = 0; s < numSupp; s++, offset++) { 828 for (q = 0; q < Nq * dE * dE; q++) { 829 geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q]; 830 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 831 } 832 for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q]; 833 } 834 } 835 PetscCall(PetscFEGeomDestroy(&cellGeom)); 836 PetscCall(ISDestroy(&suppIS)); 837 PetscCall(PetscFree(cells)); 838 } else { 839 DMField_DS *dsfield = (DMField_DS *)field->data; 840 PetscObject faceDisc, cellDisc; 841 PetscClassId faceId, cellId; 842 PetscDualSpace dsp; 843 DM K; 844 DMPolytopeType ct; 845 PetscInt(*co)[2][3]; 846 PetscInt coneSize; 847 PetscInt **counts; 848 PetscInt f, i, o, q, s; 849 const PetscInt *coneK; 850 PetscInt eStart, minOrient, maxOrient, numOrient; 851 PetscInt *orients; 852 PetscReal **orientPoints; 853 PetscReal *cellPoints; 854 PetscReal *dummyWeights; 855 PetscQuadrature cellQuad = NULL; 856 857 PetscCall(DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc)); 858 PetscCall(DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc)); 859 PetscCall(PetscObjectGetClassId(faceDisc, &faceId)); 860 PetscCall(PetscObjectGetClassId(cellDisc, &cellId)); 861 PetscCheck(faceId == PETSCFE_CLASSID && cellId == PETSCFE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported"); 862 PetscCall(PetscFEGetDualSpace((PetscFE)cellDisc, &dsp)); 863 PetscCall(PetscDualSpaceGetDM(dsp, &K)); 864 PetscCall(DMPlexGetHeightStratum(K, 1, &eStart, NULL)); 865 PetscCall(DMPlexGetCellType(K, eStart, &ct)); 866 PetscCall(DMPlexGetConeSize(K, 0, &coneSize)); 867 PetscCall(DMPlexGetCone(K, 0, &coneK)); 868 PetscCall(PetscMalloc2(numFaces, &co, coneSize, &counts)); 869 PetscCall(PetscMalloc1(dE * Nq, &cellPoints)); 870 PetscCall(PetscMalloc1(Nq, &dummyWeights)); 871 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad)); 872 PetscCall(PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights)); 873 minOrient = PETSC_MAX_INT; 874 maxOrient = PETSC_MIN_INT; 875 for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */ 876 PetscInt point = points[p]; 877 PetscInt numSupp, numChildren; 878 const PetscInt *supp; 879 880 PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, NULL)); 881 PetscCheck(!numChildren, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 882 PetscCall(DMPlexGetSupportSize(dm, point, &numSupp)); 883 PetscCheck(numSupp <= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has %" PetscInt_FMT " support, expected at most 2", point, numSupp); 884 PetscCall(DMPlexGetSupport(dm, point, &supp)); 885 for (s = 0; s < numSupp; s++) { 886 PetscInt cell = supp[s]; 887 PetscInt numCone; 888 const PetscInt *cone, *orient; 889 890 PetscCall(DMPlexGetConeSize(dm, cell, &numCone)); 891 PetscCheck(numCone == coneSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element"); 892 PetscCall(DMPlexGetCone(dm, cell, &cone)); 893 PetscCall(DMPlexGetConeOrientation(dm, cell, &orient)); 894 for (f = 0; f < coneSize; f++) { 895 if (cone[f] == point) break; 896 } 897 co[p][s][0] = f; 898 co[p][s][1] = orient[f]; 899 co[p][s][2] = cell; 900 minOrient = PetscMin(minOrient, orient[f]); 901 maxOrient = PetscMax(maxOrient, orient[f]); 902 } 903 for (; s < 2; s++) { 904 co[p][s][0] = -1; 905 co[p][s][1] = -1; 906 co[p][s][2] = -1; 907 } 908 } 909 numOrient = maxOrient + 1 - minOrient; 910 PetscCall(DMPlexGetCone(K, 0, &coneK)); 911 /* count all (face,orientation) doubles that appear */ 912 PetscCall(PetscCalloc2(numOrient, &orients, numOrient, &orientPoints)); 913 for (f = 0; f < coneSize; f++) PetscCall(PetscCalloc1(numOrient + 1, &counts[f])); 914 for (p = 0; p < numFaces; p++) { 915 for (s = 0; s < 2; s++) { 916 if (co[p][s][0] >= 0) { 917 counts[co[p][s][0]][co[p][s][1] - minOrient]++; 918 orients[co[p][s][1] - minOrient]++; 919 } 920 } 921 } 922 for (o = 0; o < numOrient; o++) { 923 if (orients[o]) { 924 PetscInt orient = o + minOrient; 925 PetscInt q; 926 927 PetscCall(PetscMalloc1(Nq * dim, &orientPoints[o])); 928 /* rotate the quadrature points appropriately */ 929 switch (ct) { 930 case DM_POLYTOPE_POINT: 931 break; 932 case DM_POLYTOPE_SEGMENT: 933 if (orient == -2 || orient == 1) { 934 for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q]; 935 } else { 936 for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q]; 937 } 938 break; 939 case DM_POLYTOPE_TRIANGLE: 940 for (q = 0; q < Nq; q++) { 941 PetscReal lambda[3]; 942 PetscReal lambdao[3]; 943 944 /* convert to barycentric */ 945 lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.; 946 lambda[1] = (geom->xi[2 * q] + 1.) / 2.; 947 lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.; 948 if (orient >= 0) { 949 for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3]; 950 } else { 951 for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3]; 952 } 953 /* convert to coordinates */ 954 orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1]; 955 orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2]; 956 } 957 break; 958 case DM_POLYTOPE_QUADRILATERAL: 959 for (q = 0; q < Nq; q++) { 960 PetscReal xi[2], xio[2]; 961 PetscInt oabs = (orient >= 0) ? orient : -(orient + 1); 962 963 xi[0] = geom->xi[2 * q]; 964 xi[1] = geom->xi[2 * q + 1]; 965 switch (oabs) { 966 case 1: 967 xio[0] = xi[1]; 968 xio[1] = -xi[0]; 969 break; 970 case 2: 971 xio[0] = -xi[0]; 972 xio[1] = -xi[1]; 973 break; 974 case 3: 975 xio[0] = -xi[1]; 976 xio[1] = xi[0]; 977 break; 978 case 0: 979 default: 980 xio[0] = xi[0]; 981 xio[1] = xi[1]; 982 break; 983 } 984 if (orient < 0) xio[0] = -xio[0]; 985 orientPoints[o][2 * q + 0] = xio[0]; 986 orientPoints[o][2 * q + 1] = xio[1]; 987 } 988 break; 989 default: 990 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]); 991 } 992 } 993 } 994 for (f = 0; f < coneSize; f++) { 995 PetscInt face = coneK[f]; 996 PetscReal v0[3]; 997 PetscReal J[9], detJ; 998 PetscInt numCells, offset; 999 PetscInt *cells; 1000 IS suppIS; 1001 1002 PetscCall(DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ)); 1003 for (o = 0; o <= numOrient; o++) { 1004 PetscFEGeom *cellGeom; 1005 1006 if (!counts[f][o]) continue; 1007 /* If this (face,orientation) double appears, 1008 * convert the face quadrature points into volume quadrature points */ 1009 for (q = 0; q < Nq; q++) { 1010 PetscReal xi0[3] = {-1., -1., -1.}; 1011 1012 CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]); 1013 } 1014 for (p = 0, numCells = 0; p < numFaces; p++) { 1015 for (s = 0; s < 2; s++) { 1016 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++; 1017 } 1018 } 1019 PetscCall(PetscMalloc1(numCells, &cells)); 1020 for (p = 0, offset = 0; p < numFaces; p++) { 1021 for (s = 0; s < 2; s++) { 1022 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2]; 1023 } 1024 } 1025 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS)); 1026 PetscCall(DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom)); 1027 for (p = 0, offset = 0; p < numFaces; p++) { 1028 for (s = 0; s < 2; s++) { 1029 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) { 1030 for (q = 0; q < Nq * dE * dE; q++) { 1031 geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q]; 1032 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 1033 } 1034 for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q]; 1035 offset++; 1036 } 1037 } 1038 } 1039 PetscCall(PetscFEGeomDestroy(&cellGeom)); 1040 PetscCall(ISDestroy(&suppIS)); 1041 PetscCall(PetscFree(cells)); 1042 } 1043 } 1044 for (o = 0; o < numOrient; o++) { 1045 if (orients[o]) PetscCall(PetscFree(orientPoints[o])); 1046 } 1047 PetscCall(PetscFree2(orients, orientPoints)); 1048 PetscCall(PetscQuadratureDestroy(&cellQuad)); 1049 for (f = 0; f < coneSize; f++) PetscCall(PetscFree(counts[f])); 1050 PetscCall(PetscFree2(co, counts)); 1051 } 1052 PetscCall(ISRestoreIndices(pointIS, &points)); 1053 PetscCall(ISDestroy(&cellIS)); 1054 PetscFunctionReturn(PETSC_SUCCESS); 1055 } 1056 1057 static PetscErrorCode DMFieldInitialize_DS(DMField field) 1058 { 1059 PetscFunctionBegin; 1060 field->ops->destroy = DMFieldDestroy_DS; 1061 field->ops->evaluate = DMFieldEvaluate_DS; 1062 field->ops->evaluateFE = DMFieldEvaluateFE_DS; 1063 field->ops->evaluateFV = DMFieldEvaluateFV_DS; 1064 field->ops->getDegree = DMFieldGetDegree_DS; 1065 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS; 1066 field->ops->view = DMFieldView_DS; 1067 field->ops->computeFaceData = DMFieldComputeFaceData_DS; 1068 PetscFunctionReturn(PETSC_SUCCESS); 1069 } 1070 1071 PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field) 1072 { 1073 DMField_DS *dsfield; 1074 1075 PetscFunctionBegin; 1076 PetscCall(PetscNew(&dsfield)); 1077 field->data = dsfield; 1078 PetscCall(DMFieldInitialize_DS(field)); 1079 PetscFunctionReturn(PETSC_SUCCESS); 1080 } 1081 1082 PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field) 1083 { 1084 DMField b; 1085 DMField_DS *dsfield; 1086 PetscObject disc = NULL, discDG = NULL; 1087 PetscSection section; 1088 PetscBool isContainer = PETSC_FALSE; 1089 PetscClassId id = -1; 1090 PetscInt numComponents = -1, dsNumFields; 1091 1092 PetscFunctionBegin; 1093 PetscCall(DMGetLocalSection(dm, §ion)); 1094 PetscCall(PetscSectionGetFieldComponents(section, fieldNum, &numComponents)); 1095 PetscCall(DMGetNumFields(dm, &dsNumFields)); 1096 if (dsNumFields) PetscCall(DMGetField(dm, fieldNum, NULL, &disc)); 1097 if (dsNumFields && dmDG) { 1098 PetscCall(DMGetField(dmDG, fieldNum, NULL, &discDG)); 1099 PetscCall(PetscObjectReference(discDG)); 1100 } 1101 if (disc) { 1102 PetscCall(PetscObjectGetClassId(disc, &id)); 1103 isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE; 1104 } 1105 if (!disc || isContainer) { 1106 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 1107 PetscFE fe; 1108 DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN; 1109 PetscInt dim, cStart, cEnd, cellHeight; 1110 1111 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 1112 PetscCall(DMGetDimension(dm, &dim)); 1113 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 1114 if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &locct)); 1115 PetscCallMPI(MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm)); 1116 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe)); 1117 PetscCall(PetscFEViewFromOptions(fe, NULL, "-field_fe_view")); 1118 disc = (PetscObject)fe; 1119 } else PetscCall(PetscObjectReference(disc)); 1120 PetscCall(PetscObjectGetClassId(disc, &id)); 1121 if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)disc, &numComponents)); 1122 else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components"); 1123 PetscCall(DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b)); 1124 PetscCall(DMFieldSetType(b, DMFIELDDS)); 1125 dsfield = (DMField_DS *)b->data; 1126 dsfield->fieldNum = fieldNum; 1127 PetscCall(DMGetDimension(dm, &dsfield->height)); 1128 dsfield->height++; 1129 PetscCall(PetscCalloc1(dsfield->height, &dsfield->disc)); 1130 dsfield->disc[0] = disc; 1131 PetscCall(PetscObjectReference((PetscObject)vec)); 1132 dsfield->vec = vec; 1133 if (dmDG) { 1134 dsfield->dmDG = dmDG; 1135 PetscCall(PetscCalloc1(dsfield->height, &dsfield->discDG)); 1136 dsfield->discDG[0] = discDG; 1137 PetscCall(PetscObjectReference((PetscObject)vecDG)); 1138 dsfield->vecDG = vecDG; 1139 } 1140 *field = b; 1141 PetscFunctionReturn(PETSC_SUCCESS); 1142 } 1143 1144 PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field) 1145 { 1146 PetscFunctionBegin; 1147 PetscCall(DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field)); 1148 PetscFunctionReturn(PETSC_SUCCESS); 1149 } 1150