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 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 nc = field->numComponents; 209 ierr = DMGetDefaultSection(field->dm,§ion);CHKERRQ(ierr); 210 ierr = DMFieldDSGetHeightDisc(field,0,&cellDisc);CHKERRQ(ierr); 211 ierr = PetscObjectGetClassId(cellDisc, &discID);CHKERRQ(ierr); 212 if (discID != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Discretization type not supported\n"); 213 cellFE = (PetscFE) cellDisc; 214 ierr = PetscFEGetDimension(cellFE,&feDim);CHKERRQ(ierr); 215 ierr = DMGetCoordinateDim(field->dm, &dim);CHKERRQ(ierr); 216 ierr = DMGetDimension(field->dm, &dimR);CHKERRQ(ierr); 217 ierr = DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF);CHKERRQ(ierr); 218 ierr = PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells);CHKERRQ(ierr); 219 for (c = 0; c < nFound; c++) { 220 if (cells[c].index < 0) SETERRQ1(PetscObjectComm((PetscObject)points),PETSC_ERR_ARG_WRONG, "Point %D could not be located\n", c); 221 } 222 ierr = PetscSFComputeDegreeBegin(cellSF,&cellDegrees);CHKERRQ(ierr); 223 ierr = PetscSFComputeDegreeEnd(cellSF,&cellDegrees);CHKERRQ(ierr); 224 for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) { 225 gatherMax = PetscMax(gatherMax,cellDegrees[c]); 226 gatherSize += cellDegrees[c]; 227 } 228 ierr = PetscMalloc3(gatherSize*dim,&cellPoints,gatherMax*dim,&coordsReal,gatherMax*dimR,&coordsRef);CHKERRQ(ierr); 229 if (datatype == PETSC_SCALAR) { 230 ierr = PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs);CHKERRQ(ierr); 231 } else { 232 ierr = PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr);CHKERRQ(ierr); 233 } 234 235 ierr = MPI_Type_contiguous(dim,MPIU_SCALAR,&pointType);CHKERRQ(ierr); 236 ierr = MPI_Type_commit(&pointType);CHKERRQ(ierr); 237 ierr = VecGetArrayRead(points,&pointsArray);CHKERRQ(ierr); 238 ierr = PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints);CHKERRQ(ierr); 239 ierr = PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints);CHKERRQ(ierr); 240 ierr = VecRestoreArrayRead(points,&pointsArray);CHKERRQ(ierr); 241 for (c = 0, offset = 0; c < numCells; c++) { 242 PetscInt nq = cellDegrees[c], p; 243 244 if (nq) { 245 PetscReal *fB, *fD, *fH; 246 PetscInt closureSize; 247 PetscScalar *elem = NULL; 248 249 for (p = 0; p < dim * nq; p++) coordsReal[p] = cellPoints[dim * offset + p]; 250 ierr = DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef);CHKERRQ(ierr); 251 ierr = PetscFEGetTabulation(cellFE,nq,coordsRef,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 252 ierr = DMPlexVecGetClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 253 if (B) { 254 if (datatype == PETSC_SCALAR) { 255 PetscScalar *cB = &cellBs[nc * offset]; 256 257 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar)); 258 } else { 259 PetscReal *cB = &cellBr[nc * offset]; 260 261 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,PetscRealPart); 262 } 263 } 264 if (D) { 265 if (datatype == PETSC_SCALAR) { 266 PetscScalar *cD = &cellDs[nc * dim * offset]; 267 268 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar)); 269 } else { 270 PetscReal *cD = &cellDr[nc * dim * offset]; 271 272 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),PetscRealPart); 273 } 274 } 275 if (H) { 276 if (datatype == PETSC_SCALAR) { 277 PetscScalar *cH = &cellHs[nc * dim * dim * offset]; 278 279 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar)); 280 } else { 281 PetscReal *cH = &cellHr[nc * dim * dim * offset]; 282 283 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart); 284 } 285 } 286 ierr = DMPlexVecRestoreClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 287 } 288 offset += nq; 289 } 290 { 291 MPI_Datatype origtype; 292 if (datatype == PETSC_SCALAR) { 293 origtype = MPIU_SCALAR; 294 } else { 295 origtype = MPIU_REAL; 296 } 297 if (B) { 298 MPI_Datatype Btype; 299 300 ierr = MPI_Type_contiguous(origtype, nc, &Btype);CHKERRQ(ierr); 301 ierr = MPI_Type_commit(&Btype);CHKERRQ(ierr); 302 ierr = PetscSFScatterBegin(cellSF,Btype,(datatype == PETSC_SCALAR) ? cellBs : cellBr, B);CHKERRQ(ierr); 303 ierr = PetscSFScatterEnd(cellSF,Btype,(datatype == PETSC_SCALAR) ? cellBs : cellBr, B);CHKERRQ(ierr); 304 ierr = MPI_Type_free(&Btype);CHKERRQ(ierr); 305 } 306 if (D) { 307 MPI_Datatype Dtype; 308 309 ierr = MPI_Type_contiguous(origtype, nc * dim, &Dtype);CHKERRQ(ierr); 310 ierr = MPI_Type_commit(&Dtype);CHKERRQ(ierr); 311 ierr = PetscSFScatterBegin(cellSF,Dtype,(datatype == PETSC_SCALAR) ? cellDs : cellDr, D);CHKERRQ(ierr); 312 ierr = PetscSFScatterEnd(cellSF,Dtype,(datatype == PETSC_SCALAR) ? cellDs : cellDr, D);CHKERRQ(ierr); 313 ierr = MPI_Type_free(&Dtype);CHKERRQ(ierr); 314 } 315 if (H) { 316 MPI_Datatype Htype; 317 318 ierr = MPI_Type_contiguous(origtype, nc * dim * dim, &Htype);CHKERRQ(ierr); 319 ierr = MPI_Type_commit(&Htype);CHKERRQ(ierr); 320 ierr = PetscSFScatterBegin(cellSF,Htype,(datatype == PETSC_SCALAR) ? cellHs : cellHr, H);CHKERRQ(ierr); 321 ierr = PetscSFScatterEnd(cellSF,Htype,(datatype == PETSC_SCALAR) ? cellHs : cellHr, H);CHKERRQ(ierr); 322 ierr = MPI_Type_free(&Htype);CHKERRQ(ierr); 323 } 324 } 325 ierr = PetscFree3(cellBr, cellDr, cellHr);CHKERRQ(ierr); 326 ierr = PetscFree3(cellBs, cellDs, cellHs);CHKERRQ(ierr); 327 ierr = PetscFree3(cellPoints,coordsReal,coordsRef);CHKERRQ(ierr); 328 ierr = MPI_Type_free(&pointType);CHKERRQ(ierr); 329 ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332 333 static PetscErrorCode DMFieldGetFEInvariance_DS(DMField field, IS pointIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 334 { 335 DMField_DS *dsfield; 336 PetscObject disc; 337 PetscInt h, imin; 338 PetscClassId id; 339 PetscErrorCode ierr; 340 341 PetscFunctionBegin; 342 dsfield = (DMField_DS *) field->data; 343 ierr = ISGetMinMax(pointIS,&imin,NULL);CHKERRQ(ierr); 344 for (h = 0; h < dsfield->height; h++) { 345 PetscInt hEnd; 346 347 ierr = DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);CHKERRQ(ierr); 348 if (imin < hEnd) break; 349 } 350 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 351 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 352 if (id == PETSCFE_CLASSID) { 353 PetscFE fe = (PetscFE) disc; 354 PetscInt order, maxOrder; 355 PetscBool tensor = PETSC_FALSE; 356 PetscSpace sp; 357 358 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 359 ierr = PetscSpaceGetOrder(sp,&order);CHKERRQ(ierr); 360 ierr = PetscSpacePolynomialGetTensor(sp,&tensor);CHKERRQ(ierr); 361 if (tensor) { 362 PetscInt dim; 363 364 ierr = DMGetDimension(field->dm,&dim);CHKERRQ(ierr); 365 maxOrder = order * dim; 366 } else { 367 maxOrder = order; 368 } 369 if (isConstant) *isConstant = (maxOrder < 1) ? PETSC_TRUE : PETSC_FALSE; 370 if (isAffine) *isAffine = (maxOrder < 2) ? PETSC_TRUE : PETSC_FALSE; 371 if (isQuadratic) *isQuadratic = (maxOrder < 3) ? PETSC_TRUE : PETSC_FALSE; 372 } 373 PetscFunctionReturn(0); 374 } 375 376 static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad) 377 { 378 PetscInt h, dim, imax, imin; 379 DM dm; 380 DMField_DS *dsfield; 381 PetscObject disc; 382 PetscFE fe; 383 PetscClassId id; 384 PetscErrorCode ierr; 385 386 387 PetscFunctionBegin; 388 dm = field->dm; 389 dsfield = (DMField_DS *) field->data; 390 ierr = ISGetMinMax(pointIS,&imax,&imin);CHKERRQ(ierr); 391 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 392 for (h = 0; h <= dim; h++) { 393 PetscInt hStart, hEnd; 394 395 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 396 if (imin >= hStart && imax < hEnd) break; 397 } 398 *quad = NULL; 399 if (h < dsfield->height) { 400 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 401 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 402 if (id != PETSCFE_CLASSID) PetscFunctionReturn(0); 403 fe = (PetscFE) disc; 404 ierr = PetscFEGetQuadrature(fe,quad);CHKERRQ(ierr); 405 ierr = PetscObjectReference((PetscObject)*quad);CHKERRQ(ierr); 406 } 407 PetscFunctionReturn(0); 408 } 409 410 static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom) 411 { 412 const PetscInt *points; 413 PetscInt p, dim, dE, numFaces, Nq; 414 PetscBool affineCells; 415 DMLabel depthLabel; 416 IS cellIS; 417 DM dm = field->dm; 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 dim = geom->dim; 422 dE = geom->dimEmbed; 423 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 424 ierr = DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);CHKERRQ(ierr); 425 ierr = DMFieldGetFEInvariance(field,cellIS,NULL,&affineCells,NULL);CHKERRQ(ierr); 426 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 427 numFaces = geom->numCells; 428 Nq = geom->numPoints; 429 if (affineCells) { 430 PetscInt numCells, offset, *cells; 431 PetscFEGeom *cellGeom; 432 IS suppIS; 433 PetscQuadrature cellQuad = NULL; 434 435 ierr = DMFieldCreateDefaultQuadrature(field,cellIS,&cellQuad);CHKERRQ(ierr); 436 for (p = 0, numCells = 0; p < numFaces; p++) { 437 PetscInt point = points[p]; 438 PetscInt numSupp, numChildren; 439 440 ierr = DMPlexGetTreeChildren(dm, point, &numChildren, NULL); CHKERRQ(ierr); 441 if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 442 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 443 if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp); 444 numCells += numSupp; 445 } 446 ierr = PetscMalloc1(numCells, &cells);CHKERRQ(ierr); 447 for (p = 0, offset = 0; p < numFaces; p++) { 448 PetscInt point = points[p]; 449 PetscInt numSupp, s; 450 const PetscInt *supp; 451 452 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 453 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 454 for (s = 0; s < numSupp; s++, offset++) { 455 cells[offset] = supp[s]; 456 } 457 } 458 ierr = ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);CHKERRQ(ierr); 459 ierr = DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);CHKERRQ(ierr); 460 for (p = 0, offset = 0; p < numFaces; p++) { 461 PetscInt point = points[p]; 462 PetscInt numSupp, s, q; 463 const PetscInt *supp; 464 465 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 466 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 467 for (s = 0; s < numSupp; s++, offset++) { 468 for (q = 0; q < Nq * dE * dE; q++) { 469 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 470 } 471 } 472 } 473 ierr = PetscFEGeomDestroy(&cellGeom);CHKERRQ(ierr); 474 ierr = ISDestroy(&suppIS);CHKERRQ(ierr); 475 ierr = PetscFree(cells);CHKERRQ(ierr); 476 ierr = PetscQuadratureDestroy(&cellQuad);CHKERRQ(ierr); 477 } else { 478 PetscObject faceDisc, cellDisc; 479 PetscClassId faceId, cellId; 480 PetscDualSpace dsp; 481 DM K; 482 PetscInt (*co)[2][3]; 483 PetscInt coneSize; 484 PetscInt **counts; 485 PetscInt f, i, o, q, s; 486 const PetscInt *coneK; 487 PetscInt minOrient, maxOrient, numOrient; 488 PetscInt *orients; 489 PetscReal **orientPoints; 490 PetscReal *cellPoints; 491 PetscReal *dummyWeights; 492 PetscQuadrature cellQuad = NULL; 493 494 ierr = DMFieldDSGetHeightDisc(field, 1, &faceDisc);CHKERRQ(ierr); 495 ierr = DMFieldDSGetHeightDisc(field, 0, &cellDisc);CHKERRQ(ierr); 496 ierr = PetscObjectGetClassId(faceDisc,&faceId);CHKERRQ(ierr); 497 ierr = PetscObjectGetClassId(cellDisc,&cellId);CHKERRQ(ierr); 498 if (faceId != PETSCFE_CLASSID || cellId != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported\n"); 499 ierr = PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);CHKERRQ(ierr); 500 ierr = PetscDualSpaceGetDM(dsp, &K); CHKERRQ(ierr); 501 ierr = DMPlexGetConeSize(K,0,&coneSize);CHKERRQ(ierr); 502 ierr = DMPlexGetCone(K,0,&coneK);CHKERRQ(ierr); 503 ierr = PetscMalloc4(numFaces,&co,dE*Nq,&cellPoints,coneSize,&counts,Nq,&dummyWeights);CHKERRQ(ierr); 504 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)field), &cellQuad);CHKERRQ(ierr); 505 ierr = PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);CHKERRQ(ierr); 506 minOrient = PETSC_MAX_INT; 507 maxOrient = PETSC_MIN_INT; 508 for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */ 509 PetscInt point = points[p]; 510 PetscInt numSupp, numChildren; 511 const PetscInt *supp; 512 513 ierr = DMPlexGetTreeChildren(dm, point, &numChildren, NULL); CHKERRQ(ierr); 514 if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children"); 515 ierr = DMPlexGetSupportSize(dm, point,&numSupp);CHKERRQ(ierr); 516 if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp); 517 ierr = DMPlexGetSupport(dm, point, &supp);CHKERRQ(ierr); 518 for (s = 0; s < numSupp; s++) { 519 PetscInt cell = supp[s]; 520 PetscInt numCone; 521 const PetscInt *cone, *orient; 522 523 ierr = DMPlexGetConeSize(dm, cell, &numCone);CHKERRQ(ierr); 524 if (numCone != coneSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element"); 525 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 526 ierr = DMPlexGetConeOrientation(dm, cell, &orient);CHKERRQ(ierr); 527 for (f = 0; f < coneSize; f++) { 528 if (cone[f] == point) break; 529 } 530 co[p][s][0] = f; 531 co[p][s][1] = orient[f]; 532 co[p][s][2] = cell; 533 minOrient = PetscMin(minOrient, orient[f]); 534 maxOrient = PetscMin(maxOrient, orient[f]); 535 } 536 for (; s < 2; s++) { 537 co[p][s][0] = -1; 538 co[p][s][1] = -1; 539 co[p][s][2] = -1; 540 } 541 } 542 numOrient = maxOrient + 1 - minOrient; 543 ierr = DMPlexGetCone(K,0,&coneK);CHKERRQ(ierr); 544 /* count all (face,orientation) doubles that appear */ 545 ierr = PetscCalloc2(numOrient,&orients,numOrient,&orientPoints);CHKERRQ(ierr); 546 for (f = 0; f < coneSize; f++) {ierr = PetscCalloc1(numOrient, &counts[f]);CHKERRQ(ierr);} 547 for (p = 0; p < numFaces; p++) { 548 for (s = 0; s < 2; s++) { 549 if (co[p][s][0] >= 0) { 550 counts[co[p][s][0]][co[p][s][1] - minOrient]++; 551 orients[co[p][s][1] - minOrient]++; 552 } 553 } 554 } 555 for (o = 0; o < numOrient; o++) { 556 if (orients[o]) { 557 PetscInt orient = o + minOrient; 558 PetscInt q; 559 560 ierr = PetscMalloc1(Nq * dim, &orientPoints[o]);CHKERRQ(ierr); 561 /* rotate the quadrature points appropriately */ 562 switch (dim) { 563 case 0: 564 break; 565 case 1: 566 if (orient == -2 || orient == 1) { 567 for (q = 0; q < Nq; q++) { 568 orientPoints[o][q] = -geom->xi[q]; 569 } 570 } else { 571 for (q = 0; q < Nq; q++) { 572 orientPoints[o][q] = geom->xi[q]; 573 } 574 } 575 break; 576 case 2: 577 switch (coneSize) { 578 case 3: 579 for (q = 0; q < Nq; q++) { 580 PetscReal lambda[3]; 581 PetscReal lambdao[3]; 582 583 /* convert to barycentric */ 584 lambda[0] = - (geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.; 585 lambda[1] = (geom->xi[2 * q] + 1.) / 2.; 586 lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.; 587 if (orient >= 0) { 588 for (i = 0; i < 3; i++) { 589 lambdao[i] = lambda[(orient + i) % 3]; 590 } 591 } else { 592 for (i = 0; i < 3; i++) { 593 lambdao[i] = lambda[(-(orient + i) + 3) % 3]; 594 } 595 } 596 /* convert to coordinates */ 597 orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1]; 598 orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2]; 599 } 600 break; 601 case 4: 602 for (q = 0; q < Nq; q++) { 603 PetscReal xi[2], xio[2]; 604 PetscInt oabs = (orient >= 0) ? orient : -(orient + 1); 605 606 xi[0] = geom->xi[2 * q]; 607 xi[1] = geom->xi[2 * q + 1]; 608 switch (oabs) { 609 case 0: 610 xio[0] = xi[0]; 611 xio[1] = xi[1]; 612 break; 613 case 1: 614 xio[0] = xi[1]; 615 xio[1] = -xi[0]; 616 break; 617 case 2: 618 xio[0] = -xi[0]; 619 xio[1] = -xi[1]; 620 case 3: 621 xio[0] = -xi[1]; 622 xio[1] = xi[0]; 623 } 624 if (orient < 0) { 625 xio[0] = -xio[0]; 626 } 627 orientPoints[o][2 * q + 0] = xio[0]; 628 orientPoints[o][2 * q + 1] = xio[1]; 629 } 630 break; 631 default: 632 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not implemented yet\n"); 633 } 634 default: 635 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not implemented yet\n"); 636 } 637 } 638 } 639 for (f = 0; f < coneSize; f++) { 640 PetscInt face = coneK[f]; 641 PetscScalar v0[3]; 642 PetscScalar J[9]; 643 PetscInt numCells, offset; 644 PetscInt *cells; 645 IS suppIS; 646 647 ierr = DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, NULL);CHKERRQ(ierr); 648 for (o = 0; o <= numOrient; o++) { 649 PetscFEGeom *cellGeom; 650 651 if (!counts[f][o]) continue; 652 /* If this (face,orientation) double appears, 653 * convert the face quadrature points into volume quadrature points */ 654 for (q = 0; q < Nq; q++) { 655 PetscReal xi0[3] = {-1., -1., -1.}; 656 657 CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]); 658 } 659 for (p = 0, numCells = 0; p < numFaces; p++) { 660 for (s = 0; s < 2; s++) { 661 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++; 662 } 663 } 664 ierr = PetscMalloc1(numCells, &cells);CHKERRQ(ierr); 665 for (p = 0, offset = 0; p < numFaces; p++) { 666 for (s = 0; s < 2; s++) { 667 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) { 668 cells[offset++] = co[p][s][2]; 669 } 670 } 671 } 672 ierr = ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);CHKERRQ(ierr); 673 ierr = DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);CHKERRQ(ierr); 674 for (p = 0, offset = 0; p < numFaces; p++) { 675 for (s = 0; s < 2; s++) { 676 if (co[p][s][0] == f && co[p][s][1] == o + minOrient) { 677 for (q = 0; q < Nq * dE * dE; q++) { 678 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q]; 679 } 680 offset++; 681 } 682 } 683 } 684 ierr = PetscFEGeomDestroy(&cellGeom);CHKERRQ(ierr); 685 ierr = ISDestroy(&suppIS);CHKERRQ(ierr); 686 ierr = PetscFree(cells);CHKERRQ(ierr); 687 } 688 } 689 for (o = 0; o < numOrient; o++) { 690 if (orients[o]) { 691 ierr = PetscFree(orientPoints[o]);CHKERRQ(ierr); 692 } 693 } 694 ierr = PetscFree2(orients,orientPoints);CHKERRQ(ierr); 695 ierr = PetscQuadratureDestroy(&cellQuad);CHKERRQ(ierr); 696 ierr = PetscFree4(co,cellPoints,counts,dummyWeights);CHKERRQ(ierr); 697 } 698 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 699 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 700 PetscFunctionReturn(0); 701 } 702 703 static PetscErrorCode DMFieldInitialize_DS(DMField field) 704 { 705 PetscFunctionBegin; 706 field->ops->destroy = DMFieldDestroy_DS; 707 field->ops->evaluate = DMFieldEvaluate_DS; 708 field->ops->evaluateFE = DMFieldEvaluateFE_DS; 709 field->ops->getFEInvariance = DMFieldGetFEInvariance_DS; 710 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS; 711 field->ops->view = DMFieldView_DS; 712 field->ops->computeFaceData = DMFieldComputeFaceData_DS; 713 PetscFunctionReturn(0); 714 } 715 716 PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field) 717 { 718 DMField_DS *dsfield; 719 PetscErrorCode ierr; 720 721 PetscFunctionBegin; 722 ierr = PetscNewLog(field,&dsfield);CHKERRQ(ierr); 723 field->data = dsfield; 724 ierr = DMFieldInitialize_DS(field);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field) 729 { 730 DMField b; 731 DMField_DS *dsfield; 732 PetscObject disc = NULL; 733 PetscBool isContainer = PETSC_FALSE; 734 PetscClassId id = -1; 735 PetscInt numComponents = -1, dsNumFields; 736 PetscSection section; 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 741 ierr = PetscSectionGetFieldComponents(section,fieldNum,&numComponents);CHKERRQ(ierr); 742 ierr = DMGetNumFields(dm,&dsNumFields);CHKERRQ(ierr); 743 if (dsNumFields) {ierr = DMGetField(dm,fieldNum,&disc);CHKERRQ(ierr);} 744 if (disc) { 745 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 746 isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE; 747 } 748 if (!disc || isContainer) { 749 MPI_Comm comm = PetscObjectComm((PetscObject) dm); 750 PetscInt cStart, cEnd, dim; 751 PetscInt localConeSize = 0, coneSize; 752 PetscFE fe; 753 PetscDualSpace Q; 754 PetscSpace P; 755 DM K; 756 PetscQuadrature quad, fquad; 757 PetscBool isSimplex; 758 759 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 760 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 761 if (cEnd > cStart) { 762 ierr = DMPlexGetConeSize(dm, cStart, &localConeSize);CHKERRQ(ierr); 763 } 764 ierr = MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 765 isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE; 766 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 767 ierr = PetscSpaceSetOrder(P, 1);CHKERRQ(ierr); 768 ierr = PetscSpaceSetNumComponents(P, numComponents);CHKERRQ(ierr); 769 ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 770 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 771 ierr = PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 772 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 773 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 774 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 775 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 776 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 777 ierr = DMDestroy(&K);CHKERRQ(ierr); 778 ierr = PetscDualSpaceSetNumComponents(Q, numComponents);CHKERRQ(ierr); 779 ierr = PetscDualSpaceSetOrder(Q, 1);CHKERRQ(ierr); 780 ierr = PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 781 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 782 ierr = PetscFECreate(comm, &fe);CHKERRQ(ierr); 783 ierr = PetscFESetType(fe,PETSCFEBASIC);CHKERRQ(ierr); 784 ierr = PetscFESetBasisSpace(fe, P);CHKERRQ(ierr); 785 ierr = PetscFESetDualSpace(fe, Q);CHKERRQ(ierr); 786 ierr = PetscFESetNumComponents(fe, numComponents);CHKERRQ(ierr); 787 ierr = PetscFESetUp(fe);CHKERRQ(ierr); 788 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 789 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 790 if (isSimplex) { 791 ierr = PetscDTGaussJacobiQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 792 ierr = PetscDTGaussJacobiQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 793 } 794 else { 795 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 796 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 797 } 798 ierr = PetscFESetQuadrature(fe, quad);CHKERRQ(ierr); 799 ierr = PetscFESetFaceQuadrature(fe, fquad);CHKERRQ(ierr); 800 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 801 ierr = PetscQuadratureDestroy(&fquad);CHKERRQ(ierr); 802 disc = (PetscObject) fe; 803 } else { 804 ierr = PetscObjectReference(disc);CHKERRQ(ierr); 805 } 806 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 807 if (id == PETSCFE_CLASSID) { 808 PetscFE fe = (PetscFE) disc; 809 810 ierr = PetscFEGetNumComponents(fe,&numComponents);CHKERRQ(ierr); 811 } else {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");} 812 ierr = DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 813 ierr = DMFieldSetType(b,DMFIELDDS);CHKERRQ(ierr); 814 dsfield = (DMField_DS *) b->data; 815 dsfield->fieldNum = fieldNum; 816 ierr = DMGetDimension(dm,&dsfield->height);CHKERRQ(ierr); 817 dsfield->height++; 818 ierr = PetscCalloc1(dsfield->height,&dsfield->disc);CHKERRQ(ierr); 819 dsfield->disc[0] = disc; 820 ierr = PetscObjectReference((PetscObject)vec);CHKERRQ(ierr); 821 dsfield->vec = vec; 822 *field = b; 823 824 PetscFunctionReturn(0); 825 } 826