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