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