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