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