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