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