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