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