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