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