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