Actual source code: dmfieldds.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  1:  #include <petsc/private/dmfieldimpl.h>
  2:  #include <petsc/private/petscfeimpl.h>
  3:  #include <petscfe.h>
  4:  #include <petscdmplex.h>
  5:  #include <petscds.h>

  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;

 17: static PetscErrorCode DMFieldDestroy_DS(DMField field)
 18: {
 19:   DMField_DS     *dsfield;
 20:   PetscInt       i;

 24:   dsfield = (DMField_DS *) field->data;
 25:   VecDestroy(&dsfield->vec);
 26:   for (i = 0; i < dsfield->height; i++) {
 27:     PetscObjectDereference(dsfield->disc[i]);
 28:   }
 29:   PetscFree(dsfield->disc);
 30:   PetscFree(dsfield);
 31:   return(0);
 32: }

 34: static PetscErrorCode DMFieldView_DS(DMField field,PetscViewer viewer)
 35: {
 36:   DMField_DS     *dsfield = (DMField_DS *) field->data;
 37:   PetscBool      iascii;
 38:   PetscObject    disc;

 42:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 43:   disc = dsfield->disc[0];
 44:   if (iascii) {
 45:     PetscViewerASCIIPrintf(viewer, "PetscDS field %D\n",dsfield->fieldNum);
 46:     PetscViewerASCIIPushTab(viewer);
 47:     PetscObjectView(disc,viewer);
 48:     PetscViewerASCIIPopTab(viewer);
 49:   }
 50:   PetscViewerASCIIPushTab(viewer);
 51:   if (dsfield->multifieldVec) {
 52:     SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"View of subfield not implemented yet");
 53:   } else {
 54:     VecView(dsfield->vec,viewer);
 55:   }
 56:   PetscViewerASCIIPopTab(viewer);
 57:   return(0);
 58: }

 60: static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject *disc)
 61: {
 62:   DMField_DS     *dsfield = (DMField_DS *) field->data;

 66:   if (!dsfield->disc[height]) {
 67:     PetscClassId   id;

 69:     PetscObjectGetClassId(dsfield->disc[0],&id);
 70:     if (id == PETSCFE_CLASSID) {
 71:       PetscFE fe = (PetscFE) dsfield->disc[0];

 73:       PetscFECreateHeightTrace(fe,height,(PetscFE *)&dsfield->disc[height]);
 74:     }
 75:   }
 76:   *disc = dsfield->disc[height];
 77:   return(0);
 78: }

 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)

 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;

110:   dm   = field->dm;
111:   nc   = field->numComponents;
112:   PetscQuadratureGetData(quad,&dim,NULL,&nq,&qpoints,NULL);
113:   DMFieldDSGetHeightDisc(field,dsfield->height - 1 - dim,&disc);
114:   DMGetDimension(dm,&meshDim);
115:   DMGetLocalSection(dm,&section);
116:   PetscSectionGetField(section,dsfield->fieldNum,&section);
117:   PetscObjectGetClassId(disc,&classid);
118:   /* TODO: batch */
119:   PetscObjectTypeCompare((PetscObject)pointIS,ISSTRIDE,&isStride);
120:   ISGetLocalSize(pointIS,&numCells);
121:   if (isStride) {
122:     ISStrideGetInfo(pointIS,&sfirst,&stride);
123:   } else {
124:     ISGetIndices(pointIS,&points);
125:   }
126:   if (classid == PETSCFE_CLASSID) {
127:     PetscFE      fe = (PetscFE) disc;
128:     PetscInt     feDim, i;
129:     PetscReal    *fB = NULL, *fD = NULL, *fH = NULL;

131:     PetscFEGetDimension(fe,&feDim);
132:     PetscFEGetTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);
133:     for (i = 0; i < numCells; i++) {
134:       PetscInt     c = isStride ? (sfirst + i * stride) : points[i];
135:       PetscInt     closureSize;
136:       PetscScalar *elem = NULL;

138:       DMPlexVecGetClosure(dm,section,dsfield->vec,c,&closureSize,&elem);
139:       if (B) {
140:         if (type == PETSC_SCALAR) {
141:           PetscScalar *cB = &((PetscScalar *) B)[nc * nq * i];

143:           DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar));
144:         } else {
145:           PetscReal *cB = &((PetscReal *) B)[nc * nq * i];

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];

154:           DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar));
155:         } else {
156:           PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i];

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];

165:           DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar));
166:         } else {
167:           PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i];

169:           DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart);
170:         }
171:       }
172:       DMPlexVecRestoreClosure(dm,section,dsfield->vec,c,&closureSize,&elem);
173:     }
174:     PetscFERestoreTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);
175:   } else {SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");}
176:   if (!isStride) {
177:     ISRestoreIndices(pointIS,&points);
178:   }
179:   return(0);
180: }

182: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
183: {
184:   DMField_DS        *dsfield = (DMField_DS *) field->data;
185:   PetscSF            cellSF = NULL;
186:   const PetscSFNode *cells;
187:   PetscInt           c, nFound, numCells, feDim, nc;
188:   const PetscInt    *cellDegrees;
189:   const PetscScalar *pointsArray;
190:   PetscScalar       *cellPoints;
191:   PetscInt           gatherSize, gatherMax;
192:   PetscInt           dim, dimR, offset;
193:   MPI_Datatype       pointType;
194:   PetscObject        cellDisc;
195:   PetscFE            cellFE;
196:   PetscClassId       discID;
197:   PetscReal         *coordsReal, *coordsRef;
198:   PetscSection       section;
199:   PetscScalar       *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
200:   PetscReal         *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
201:   PetscReal         *v, *J, *invJ, *detJ;
202:   PetscErrorCode     ierr;

205:   nc   = field->numComponents;
206:   DMGetLocalSection(field->dm,&section);
207:   DMFieldDSGetHeightDisc(field,0,&cellDisc);
208:   PetscObjectGetClassId(cellDisc, &discID);
209:   if (discID != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Discretization type not supported\n");
210:   cellFE = (PetscFE) cellDisc;
211:   PetscFEGetDimension(cellFE,&feDim);
212:   DMGetCoordinateDim(field->dm, &dim);
213:   DMGetDimension(field->dm, &dimR);
214:   DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF);
215:   PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells);
216:   for (c = 0; c < nFound; c++) {
217:     if (cells[c].index < 0) SETERRQ1(PetscObjectComm((PetscObject)points),PETSC_ERR_ARG_WRONG, "Point %D could not be located\n", c);
218:   }
219:   PetscSFComputeDegreeBegin(cellSF,&cellDegrees);
220:   PetscSFComputeDegreeEnd(cellSF,&cellDegrees);
221:   for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
222:     gatherMax = PetscMax(gatherMax,cellDegrees[c]);
223:     gatherSize += cellDegrees[c];
224:   }
225:   PetscMalloc3(gatherSize*dim,&cellPoints,gatherMax*dim,&coordsReal,gatherMax*dimR,&coordsRef);
226:   PetscMalloc4(gatherMax*dimR,&v,gatherMax*dimR*dimR,&J,gatherMax*dimR*dimR,&invJ,gatherMax,&detJ);
227:   if (datatype == PETSC_SCALAR) {
228:     PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs);
229:   } else {
230:     PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr);
231:   }

233:   MPI_Type_contiguous(dim,MPIU_SCALAR,&pointType);
234:   MPI_Type_commit(&pointType);
235:   VecGetArrayRead(points,&pointsArray);
236:   PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints);
237:   PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints);
238:   VecRestoreArrayRead(points,&pointsArray);
239:   for (c = 0, offset = 0; c < numCells; c++) {
240:     PetscInt nq = cellDegrees[c], p;

242:     if (nq) {
243:       PetscReal *fB, *fD, *fH;
244:       PetscInt     closureSize;
245:       PetscScalar *elem = NULL;
246:       PetscReal   *quadPoints;
247:       PetscQuadrature quad;
248:       PetscInt d, e, f, g;

250:       for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
251:       DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef);
252:       PetscFEGetTabulation(cellFE,nq,coordsRef,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);
253:       PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
254:       PetscMalloc1(dimR * nq, &quadPoints);
255:       for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
256:       PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL);
257:       DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ);
258:       PetscQuadratureDestroy(&quad);
259:       DMPlexVecGetClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);
260:       if (B) {
261:         if (datatype == PETSC_SCALAR) {
262:           PetscScalar *cB = &cellBs[nc * offset];

264:           DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar));
265:         } else {
266:           PetscReal *cB = &cellBr[nc * offset];

268:           DMFieldDSdot(cB,fB,elem,nq,feDim,nc,PetscRealPart);
269:         }
270:       }
271:       if (D) {
272:         if (datatype == PETSC_SCALAR) {
273:           PetscScalar *cD = &cellDs[nc * dim * offset];

275:           DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar));
276:           for (p = 0; p < nq; p++) {
277:             for (g = 0; g < nc; g++) {
278:               PetscScalar vs[3];

280:               for (d = 0; d < dimR; d++) {
281:                 vs[d] = 0.;
282:                 for (e = 0; e < dimR; e++) {
283:                   vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
284:                 }
285:               }
286:               for (d = 0; d < dimR; d++) {
287:                 cD[(nc * p + g) * dimR + d] = vs[d];
288:               }
289:             }
290:           }
291:         } else {
292:           PetscReal *cD = &cellDr[nc * dim * offset];

294:           DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),PetscRealPart);
295:           for (p = 0; p < nq; p++) {
296:             for (g = 0; g < nc; g++) {
297:               for (d = 0; d < dimR; d++) {
298:                 v[d] = 0.;
299:                 for (e = 0; e < dimR; e++) {
300:                   v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
301:                 }
302:               }
303:               for (d = 0; d < dimR; d++) {
304:                 cD[(nc * p + g) * dimR + d] = v[d];
305:               }
306:             }
307:           }
308:         }
309:       }
310:       if (H) {
311:         if (datatype == PETSC_SCALAR) {
312:           PetscScalar *cH = &cellHs[nc * dim * dim * offset];

314:           DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar));
315:           for (p = 0; p < nq; p++) {
316:             for (g = 0; g < nc * dimR; g++) {
317:               PetscScalar vs[3];

319:               for (d = 0; d < dimR; d++) {
320:                 vs[d] = 0.;
321:                 for (e = 0; e < dimR; e++) {
322:                   vs[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] = vs[d];
327:               }
328:             }
329:             for (g = 0; g < nc; g++) {
330:               for (f = 0; f < dimR; f++) {
331:                 PetscScalar vs[3];

333:                 for (d = 0; d < dimR; d++) {
334:                   vs[d] = 0.;
335:                   for (e = 0; e < dimR; e++) {
336:                     vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
337:                   }
338:                 }
339:                 for (d = 0; d < dimR; d++) {
340:                   cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
341:                 }
342:               }
343:             }
344:           }
345:         } else {
346:           PetscReal *cH = &cellHr[nc * dim * dim * offset];

348:           DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart);
349:           for (p = 0; p < nq; p++) {
350:             for (g = 0; g < nc * dimR; g++) {
351:               for (d = 0; d < dimR; d++) {
352:                 v[d] = 0.;
353:                 for (e = 0; e < dimR; e++) {
354:                   v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
355:                 }
356:               }
357:               for (d = 0; d < dimR; d++) {
358:                 cH[(nc * dimR * p + g) * dimR + d] = v[d];
359:               }
360:             }
361:             for (g = 0; g < nc; g++) {
362:               for (f = 0; f < dimR; f++) {
363:                 for (d = 0; d < dimR; d++) {
364:                   v[d] = 0.;
365:                   for (e = 0; e < dimR; e++) {
366:                     v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
367:                   }
368:                 }
369:                 for (d = 0; d < dimR; d++) {
370:                   cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
371:                 }
372:               }
373:             }
374:           }
375:         }
376:       }
377:       DMPlexVecRestoreClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);
378:       PetscFERestoreTabulation(cellFE,nq,coordsRef,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);
379:     }
380:     offset += nq;
381:   }
382:   {
383:     MPI_Datatype origtype;

385:     if (datatype == PETSC_SCALAR) {
386:       origtype = MPIU_SCALAR;
387:     } else {
388:       origtype = MPIU_REAL;
389:     }
390:     if (B) {
391:       MPI_Datatype Btype;

393:       MPI_Type_contiguous(nc, origtype, &Btype);
394:       MPI_Type_commit(&Btype);
395:       PetscSFScatterBegin(cellSF,Btype,(datatype == PETSC_SCALAR) ? (void *) cellBs : (void *) cellBr, B);
396:       PetscSFScatterEnd(cellSF,Btype,(datatype == PETSC_SCALAR) ? (void *) cellBs : (void *) cellBr, B);
397:       MPI_Type_free(&Btype);
398:     }
399:     if (D) {
400:       MPI_Datatype Dtype;

402:       MPI_Type_contiguous(nc * dim, origtype, &Dtype);
403:       MPI_Type_commit(&Dtype);
404:       PetscSFScatterBegin(cellSF,Dtype,(datatype == PETSC_SCALAR) ? (void *) cellDs : (void *) cellDr, D);
405:       PetscSFScatterEnd(cellSF,Dtype,(datatype == PETSC_SCALAR) ? (void *) cellDs : (void *) cellDr, D);
406:       MPI_Type_free(&Dtype);
407:     }
408:     if (H) {
409:       MPI_Datatype Htype;

411:       MPI_Type_contiguous(nc * dim * dim, origtype, &Htype);
412:       MPI_Type_commit(&Htype);
413:       PetscSFScatterBegin(cellSF,Htype,(datatype == PETSC_SCALAR) ? (void *) cellHs : (void *) cellHr, H);
414:       PetscSFScatterEnd(cellSF,Htype,(datatype == PETSC_SCALAR) ? (void *) cellHs : (void *) cellHr, H);
415:       MPI_Type_free(&Htype);
416:     }
417:   }
418:   PetscFree4(v,J,invJ,detJ);
419:   PetscFree3(cellBr, cellDr, cellHr);
420:   PetscFree3(cellBs, cellDs, cellHs);
421:   PetscFree3(cellPoints,coordsReal,coordsRef);
422:   MPI_Type_free(&pointType);
423:   PetscSFDestroy(&cellSF);
424:   return(0);
425: }

427: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
428: {
429:   DMField_DS      *dsfield = (DMField_DS *) field->data;
430:   PetscInt         h, imin;
431:   PetscInt         dim;
432:   PetscClassId     id;
433:   PetscQuadrature  quad = NULL;
434:   PetscInt         maxDegree;
435:   PetscFEGeom      *geom;
436:   PetscInt         Nq, Nc, dimC, qNc, N;
437:   PetscInt         numPoints;
438:   void            *qB = NULL, *qD = NULL, *qH = NULL;
439:   const PetscReal *weights;
440:   MPI_Datatype     mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
441:   PetscObject      disc;
442:   DMField          coordField;
443:   PetscErrorCode   ierr;

446:   Nc = field->numComponents;
447:   DMGetCoordinateDim(field->dm, &dimC);
448:   DMGetDimension(field->dm, &dim);
449:   ISGetLocalSize(pointIS, &numPoints);
450:   ISGetMinMax(pointIS,&imin,NULL);
451:   for (h = 0; h < dsfield->height; h++) {
452:     PetscInt hEnd;

454:     DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);
455:     if (imin < hEnd) break;
456:   }
457:   dim -= h;
458:   DMFieldDSGetHeightDisc(field,h,&disc);
459:   PetscObjectGetClassId(disc,&id);
460:   if (id != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported\n");
461:   DMGetCoordinateField(field->dm, &coordField);
462:   DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
463:   if (maxDegree <= 1) {
464:     DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
465:   }
466:   if (!quad) {DMFieldCreateDefaultQuadrature(field, pointIS, &quad);}
467:   if (!quad) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages\n");
468:   DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&geom);
469:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights);
470:   if (qNc != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components\n");
471:   N = numPoints * Nq * Nc;
472:   if (B) DMGetWorkArray(field->dm, N, mpitype, &qB);
473:   if (D) DMGetWorkArray(field->dm, N * dimC, mpitype, &qD);
474:   if (H) DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
475:   DMFieldEvaluateFE(field,pointIS,quad,type,qB,qD,qH);
476:   if (B) {
477:     PetscInt i, j, k;

479:     if (type == PETSC_SCALAR) {
480:       PetscScalar * sB  = (PetscScalar *) B;
481:       PetscScalar * sqB = (PetscScalar *) qB;

483:       for (i = 0; i < numPoints; i++) {
484:         PetscReal vol = 0.;

486:         for (j = 0; j < Nc; j++) {sB[i * Nc + j] = 0.;}
487:         for (k = 0; k < Nq; k++) {
488:           vol += geom->detJ[i * Nq + k] * weights[k];
489:           for (j = 0; j < Nc; j++) {
490:             sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[ (i * Nq + k) * Nc + j];
491:           }
492:         }
493:         for (k = 0; k < Nq * Nc; k++) sB[i * Nq * Nc + k] /= vol;
494:       }
495:     } else {
496:       PetscReal * rB  = (PetscReal *) B;
497:       PetscReal * rqB = (PetscReal *) qB;

499:       for (i = 0; i < numPoints; i++) {
500:         PetscReal vol = 0.;

502:         for (j = 0; j < Nc; j++) {rB[i * Nc + j] = 0.;}
503:         for (k = 0; k < Nq; k++) {
504:           vol += geom->detJ[i * Nq + k] * weights[k];
505:           for (j = 0; j < Nc; j++) {
506:             rB[i * Nc + j] += weights[k] * rqB[ (i * Nq + k) * Nc + j];
507:           }
508:         }
509:         for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
510:       }
511:     }
512:   }
513:   if (D) {
514:     PetscInt i, j, k, l, m;

516:     if (type == PETSC_SCALAR) {
517:       PetscScalar * sD  = (PetscScalar *) D;
518:       PetscScalar * sqD = (PetscScalar *) qD;

520:       for (i = 0; i < numPoints; i++) {
521:         PetscReal vol = 0.;

523:         for (j = 0; j < Nc * dimC; j++) {sD[i * Nc * dimC + j] = 0.;}
524:         for (k = 0; k < Nq; k++) {
525:           vol += geom->detJ[i * Nq + k] * weights[k];
526:           for (j = 0; j < Nc; j++) {
527:             PetscScalar pD[3] = {0.,0.,0.};

529:             for (l = 0; l < dimC; l++) {
530:               for (m = 0; m < dim; m++) {
531:                 pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m];
532:               }
533:             }
534:             for (l = 0; l < dimC; l++) {
535:               sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
536:             }
537:           }
538:         }
539:         for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
540:       }
541:     } else {
542:       PetscReal * rD  = (PetscReal *) D;
543:       PetscReal * rqD = (PetscReal *) qD;

545:       for (i = 0; i < numPoints; i++) {
546:         PetscReal vol = 0.;

548:         for (j = 0; j < Nc * dimC; j++) {rD[i * Nc * dimC + j] = 0.;}
549:         for (k = 0; k < Nq; k++) {
550:           vol += geom->detJ[i * Nq + k] * weights[k];
551:           for (j = 0; j < Nc; j++) {
552:             PetscReal pD[3] = {0.,0.,0.};

554:             for (l = 0; l < dimC; l++) {
555:               for (m = 0; m < dim; m++) {
556:                 pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m];
557:               }
558:             }
559:             for (l = 0; l < dimC; l++) {
560:               rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
561:             }
562:           }
563:         }
564:         for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
565:       }
566:     }
567:   }
568:   if (H) {
569:     PetscInt i, j, k, l, m, q, r;

571:     if (type == PETSC_SCALAR) {
572:       PetscScalar * sH  = (PetscScalar *) H;
573:       PetscScalar * sqH = (PetscScalar *) qH;

575:       for (i = 0; i < numPoints; i++) {
576:         PetscReal vol = 0.;

578:         for (j = 0; j < Nc * dimC * dimC; j++) {sH[i * Nc * dimC * dimC + j] = 0.;}
579:         for (k = 0; k < Nq; k++) {
580:           const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];

582:           vol += geom->detJ[i * Nq + k] * weights[k];
583:           for (j = 0; j < Nc; j++) {
584:             PetscScalar pH[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
585:             const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];

587:             for (l = 0; l < dimC; l++) {
588:               for (m = 0; m < dimC; m++) {
589:                 for (q = 0; q < dim; q++) {
590:                   for (r = 0; r < dim; r++) {
591:                     pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
592:                   }
593:                 }
594:               }
595:             }
596:             for (l = 0; l < dimC; l++) {
597:               for (m = 0; m < dimC; m++) {
598:                 sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
599:               }
600:             }
601:           }
602:         }
603:         for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
604:       }
605:     } else {
606:       PetscReal * rH  = (PetscReal *) H;
607:       PetscReal * rqH = (PetscReal *) qH;

609:       for (i = 0; i < numPoints; i++) {
610:         PetscReal vol = 0.;

612:         for (j = 0; j < Nc * dimC * dimC; j++) {rH[i * Nc * dimC * dimC + j] = 0.;}
613:         for (k = 0; k < Nq; k++) {
614:           const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];

616:           vol += geom->detJ[i * Nq + k] * weights[k];
617:           for (j = 0; j < Nc; j++) {
618:             PetscReal pH[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
619:             const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];

621:             for (l = 0; l < dimC; l++) {
622:               for (m = 0; m < dimC; m++) {
623:                 for (q = 0; q < dim; q++) {
624:                   for (r = 0; r < dim; r++) {
625:                     pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
626:                   }
627:                 }
628:               }
629:             }
630:             for (l = 0; l < dimC; l++) {
631:               for (m = 0; m < dimC; m++) {
632:                 rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
633:               }
634:             }
635:           }
636:         }
637:         for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
638:       }
639:     }
640:   }
641:   if (B) DMRestoreWorkArray(field->dm, N, mpitype, &qB);
642:   if (D) DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD);
643:   if (H) DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
644:   PetscFEGeomDestroy(&geom);
645:   PetscQuadratureDestroy(&quad);
646:   return(0);
647: }

649: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
650: {
651:   DMField_DS     *dsfield;
652:   PetscObject    disc;
653:   PetscInt       h, imin, imax;
654:   PetscClassId   id;

658:   dsfield = (DMField_DS *) field->data;
659:   ISGetMinMax(pointIS,&imin,&imax);
660:   if (imin >= imax) {
661:     h = 0;
662:   } else {
663:     for (h = 0; h < dsfield->height; h++) {
664:       PetscInt hEnd;

666:       DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);
667:       if (imin < hEnd) break;
668:     }
669:   }
670:   DMFieldDSGetHeightDisc(field,h,&disc);
671:   PetscObjectGetClassId(disc,&id);
672:   if (id == PETSCFE_CLASSID) {
673:     PetscFE    fe = (PetscFE) disc;
674:     PetscSpace sp;

676:     PetscFEGetBasisSpace(fe, &sp);
677:     PetscSpaceGetDegree(sp, minDegree, maxDegree);
678:   }
679:   return(0);
680: }

682: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
683: {
684:   PetscInt       h, dim, imax, imin, cellHeight;
685:   DM             dm;
686:   DMField_DS     *dsfield;
687:   PetscObject    disc;
688:   PetscFE        fe;
689:   PetscClassId   id;


694:   dm = field->dm;
695:   dsfield = (DMField_DS *) field->data;
696:   ISGetMinMax(pointIS,&imin,&imax);
697:   DMGetDimension(dm,&dim);
698:   for (h = 0; h <= dim; h++) {
699:     PetscInt hStart, hEnd;

701:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
702:     if (imax >= hStart && imin < hEnd) break;
703:   }
704:   DMPlexGetVTKCellHeight(dm, &cellHeight);
705:   h -= cellHeight;
706:   *quad = NULL;
707:   if (h < dsfield->height) {
708:     DMFieldDSGetHeightDisc(field,h,&disc);
709:     PetscObjectGetClassId(disc,&id);
710:     if (id != PETSCFE_CLASSID) return(0);
711:     fe = (PetscFE) disc;
712:     PetscFEGetQuadrature(fe,quad);
713:     PetscObjectReference((PetscObject)*quad);
714:   }
715:   return(0);
716: }

718: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
719: {
720:   const PetscInt *points;
721:   PetscInt        p, dim, dE, numFaces, Nq;
722:   PetscInt        maxDegree;
723:   DMLabel         depthLabel;
724:   IS              cellIS;
725:   DM              dm = field->dm;
726:   PetscErrorCode  ierr;

729:   dim = geom->dim;
730:   dE  = geom->dimEmbed;
731:   DMPlexGetDepthLabel(dm, &depthLabel);
732:   DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);
733:   DMFieldGetDegree(field,cellIS,NULL,&maxDegree);
734:   ISGetIndices(pointIS, &points);
735:   numFaces = geom->numCells;
736:   Nq = geom->numPoints;
737:   if (maxDegree <= 1) {
738:     PetscInt        numCells, offset, *cells;
739:     PetscFEGeom     *cellGeom;
740:     IS              suppIS;
741:     PetscQuadrature cellQuad = NULL;

743:     DMFieldCreateDefaultQuadrature(field,cellIS,&cellQuad);
744:     for (p = 0, numCells = 0; p < numFaces; p++) {
745:       PetscInt        point = points[p];
746:       PetscInt        numSupp, numChildren;

748:       DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
749:       if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
750:       DMPlexGetSupportSize(dm, point,&numSupp);
751:       if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp);
752:       numCells += numSupp;
753:     }
754:     PetscMalloc1(numCells, &cells);
755:     for (p = 0, offset = 0; p < numFaces; p++) {
756:       PetscInt        point = points[p];
757:       PetscInt        numSupp, s;
758:       const PetscInt *supp;

760:       DMPlexGetSupportSize(dm, point,&numSupp);
761:       DMPlexGetSupport(dm, point, &supp);
762:       for (s = 0; s < numSupp; s++, offset++) {
763:         cells[offset] = supp[s];
764:       }
765:     }
766:     ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);
767:     DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);
768:     for (p = 0, offset = 0; p < numFaces; p++) {
769:       PetscInt        point = points[p];
770:       PetscInt        numSupp, s, q;
771:       const PetscInt *supp;

773:       DMPlexGetSupportSize(dm, point,&numSupp);
774:       DMPlexGetSupport(dm, point, &supp);
775:       for (s = 0; s < numSupp; s++, offset++) {
776:         for (q = 0; q < Nq * dE * dE; q++) {
777:           geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
778:           geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
779:         }
780:         for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
781:       }
782:     }
783:     PetscFEGeomDestroy(&cellGeom);
784:     ISDestroy(&suppIS);
785:     PetscFree(cells);
786:     PetscQuadratureDestroy(&cellQuad);
787:   } else {
788:     PetscObject          faceDisc, cellDisc;
789:     PetscClassId         faceId, cellId;
790:     PetscDualSpace       dsp;
791:     DM                   K;
792:     PetscInt           (*co)[2][3];
793:     PetscInt             coneSize;
794:     PetscInt           **counts;
795:     PetscInt             f, i, o, q, s;
796:     const PetscInt      *coneK;
797:     PetscInt             minOrient, maxOrient, numOrient;
798:     PetscInt            *orients;
799:     PetscReal          **orientPoints;
800:     PetscReal           *cellPoints;
801:     PetscReal           *dummyWeights;
802:     PetscQuadrature      cellQuad = NULL;

804:     DMFieldDSGetHeightDisc(field, 1, &faceDisc);
805:     DMFieldDSGetHeightDisc(field, 0, &cellDisc);
806:     PetscObjectGetClassId(faceDisc,&faceId);
807:     PetscObjectGetClassId(cellDisc,&cellId);
808:     if (faceId != PETSCFE_CLASSID || cellId != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported\n");
809:     PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);
810:     PetscDualSpaceGetDM(dsp, &K);
811:     DMPlexGetConeSize(K,0,&coneSize);
812:     DMPlexGetCone(K,0,&coneK);
813:     PetscMalloc2(numFaces, &co, coneSize, &counts);
814:     PetscMalloc1(dE*Nq, &cellPoints);
815:     PetscMalloc1(Nq, &dummyWeights);
816:     PetscQuadratureCreate(PetscObjectComm((PetscObject)field), &cellQuad);
817:     PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);
818:     minOrient = PETSC_MAX_INT;
819:     maxOrient = PETSC_MIN_INT;
820:     for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
821:       PetscInt        point = points[p];
822:       PetscInt        numSupp, numChildren;
823:       const PetscInt *supp;

825:       DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
826:       if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
827:       DMPlexGetSupportSize(dm, point,&numSupp);
828:       if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp);
829:       DMPlexGetSupport(dm, point, &supp);
830:       for (s = 0; s < numSupp; s++) {
831:         PetscInt        cell = supp[s];
832:         PetscInt        numCone;
833:         const PetscInt *cone, *orient;

835:         DMPlexGetConeSize(dm, cell, &numCone);
836:         if (numCone != coneSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
837:         DMPlexGetCone(dm, cell, &cone);
838:         DMPlexGetConeOrientation(dm, cell, &orient);
839:         for (f = 0; f < coneSize; f++) {
840:           if (cone[f] == point) break;
841:         }
842:         co[p][s][0] = f;
843:         co[p][s][1] = orient[f];
844:         co[p][s][2] = cell;
845:         minOrient = PetscMin(minOrient, orient[f]);
846:         maxOrient = PetscMax(maxOrient, orient[f]);
847:       }
848:       for (; s < 2; s++) {
849:         co[p][s][0] = -1;
850:         co[p][s][1] = -1;
851:         co[p][s][2] = -1;
852:       }
853:     }
854:     numOrient = maxOrient + 1 - minOrient;
855:     DMPlexGetCone(K,0,&coneK);
856:     /* count all (face,orientation) doubles that appear */
857:     PetscCalloc2(numOrient,&orients,numOrient,&orientPoints);
858:     for (f = 0; f < coneSize; f++) {PetscCalloc1(numOrient+1, &counts[f]);}
859:     for (p = 0; p < numFaces; p++) {
860:       for (s = 0; s < 2; s++) {
861:         if (co[p][s][0] >= 0) {
862:           counts[co[p][s][0]][co[p][s][1] - minOrient]++;
863:           orients[co[p][s][1] - minOrient]++;
864:         }
865:       }
866:     }
867:     for (o = 0; o < numOrient; o++) {
868:       if (orients[o]) {
869:         PetscInt orient = o + minOrient;
870:         PetscInt q;

872:         PetscMalloc1(Nq * dim, &orientPoints[o]);
873:         /* rotate the quadrature points appropriately */
874:         switch (dim) {
875:         case 0:
876:           break;
877:         case 1:
878:           if (orient == -2 || orient == 1) {
879:             for (q = 0; q < Nq; q++) {
880:               orientPoints[o][q] = -geom->xi[q];
881:             }
882:           } else {
883:             for (q = 0; q < Nq; q++) {
884:               orientPoints[o][q] = geom->xi[q];
885:             }
886:           }
887:           break;
888:         case 2:
889:           switch (coneSize) {
890:           case 3:
891:             for (q = 0; q < Nq; q++) {
892:               PetscReal lambda[3];
893:               PetscReal lambdao[3];

895:               /* convert to barycentric */
896:               lambda[0] = - (geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
897:               lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
898:               lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
899:               if (orient >= 0) {
900:                 for (i = 0; i < 3; i++) {
901:                   lambdao[i] = lambda[(orient + i) % 3];
902:                 }
903:               } else {
904:                 for (i = 0; i < 3; i++) {
905:                   lambdao[i] = lambda[(-(orient + i) + 3) % 3];
906:                 }
907:               }
908:               /* convert to coordinates */
909:               orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
910:               orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
911:             }
912:             break;
913:           case 4:
914:             for (q = 0; q < Nq; q++) {
915:               PetscReal xi[2], xio[2];
916:               PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);

918:               xi[0] = geom->xi[2 * q];
919:               xi[1] = geom->xi[2 * q + 1];
920:               switch (oabs) {
921:               case 1:
922:                 xio[0] = xi[1];
923:                 xio[1] = -xi[0];
924:                 break;
925:               case 2:
926:                 xio[0] = -xi[0];
927:                 xio[1] = -xi[1];
928:               case 3:
929:                 xio[0] = -xi[1];
930:                 xio[1] = xi[0];
931:               case 0:
932:               default:
933:                 xio[0] = xi[0];
934:                 xio[1] = xi[1];
935:                 break;
936:               }
937:               if (orient < 0) {
938:                 xio[0] = -xio[0];
939:               }
940:               orientPoints[o][2 * q + 0] = xio[0];
941:               orientPoints[o][2 * q + 1] = xio[1];
942:             }
943:             break;
944:           default:
945:             SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cone size %D not yet supported\n", coneSize);
946:           }
947:         default:
948:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D not yet supported\n", dim);
949:         }
950:       }
951:     }
952:     for (f = 0; f < coneSize; f++) {
953:       PetscInt face = coneK[f];
954:       PetscReal v0[3];
955:       PetscReal J[9], detJ;
956:       PetscInt numCells, offset;
957:       PetscInt *cells;
958:       IS suppIS;

960:       DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ);
961:       for (o = 0; o <= numOrient; o++) {
962:         PetscFEGeom *cellGeom;

964:         if (!counts[f][o]) continue;
965:         /* If this (face,orientation) double appears,
966:          * convert the face quadrature points into volume quadrature points */
967:         for (q = 0; q < Nq; q++) {
968:           PetscReal xi0[3] = {-1., -1., -1.};

970:           CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
971:         }
972:         for (p = 0, numCells = 0; p < numFaces; p++) {
973:           for (s = 0; s < 2; s++) {
974:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
975:           }
976:         }
977:         PetscMalloc1(numCells, &cells);
978:         for (p = 0, offset = 0; p < numFaces; p++) {
979:           for (s = 0; s < 2; s++) {
980:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
981:               cells[offset++] = co[p][s][2];
982:             }
983:           }
984:         }
985:         ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);
986:         DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);
987:         for (p = 0, offset = 0; p < numFaces; p++) {
988:           for (s = 0; s < 2; s++) {
989:             if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
990:               for (q = 0; q < Nq * dE * dE; q++) {
991:                 geom->suppJ[s][p * Nq * dE * dE + q]    = cellGeom->J[offset * Nq * dE * dE + q];
992:                 geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
993:               }
994:               for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
995:               offset++;
996:             }
997:           }
998:         }
999:         PetscFEGeomDestroy(&cellGeom);
1000:         ISDestroy(&suppIS);
1001:         PetscFree(cells);
1002:       }
1003:     }
1004:     for (o = 0; o < numOrient; o++) {
1005:       if (orients[o]) {
1006:         PetscFree(orientPoints[o]);
1007:       }
1008:     }
1009:     PetscFree2(orients,orientPoints);
1010:     PetscQuadratureDestroy(&cellQuad);
1011:     for (f = 0; f < coneSize; f++) {PetscFree(counts[f]);}
1012:     PetscFree2(co,counts);
1013:   }
1014:   ISRestoreIndices(pointIS, &points);
1015:   ISDestroy(&cellIS);
1016:   return(0);
1017: }

1019: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1020: {
1022:   field->ops->destroy                 = DMFieldDestroy_DS;
1023:   field->ops->evaluate                = DMFieldEvaluate_DS;
1024:   field->ops->evaluateFE              = DMFieldEvaluateFE_DS;
1025:   field->ops->evaluateFV              = DMFieldEvaluateFV_DS;
1026:   field->ops->getDegree               = DMFieldGetDegree_DS;
1027:   field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1028:   field->ops->view                    = DMFieldView_DS;
1029:   field->ops->computeFaceData         = DMFieldComputeFaceData_DS;
1030:   return(0);
1031: }

1033: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1034: {
1035:   DMField_DS     *dsfield;

1039:   PetscNewLog(field,&dsfield);
1040:   field->data = dsfield;
1041:   DMFieldInitialize_DS(field);
1042:   return(0);
1043: }

1045: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field)
1046: {
1047:   DMField        b;
1048:   DMField_DS     *dsfield;
1049:   PetscObject    disc = NULL;
1050:   PetscBool      isContainer = PETSC_FALSE;
1051:   PetscClassId   id = -1;
1052:   PetscInt       numComponents = -1, dsNumFields;
1053:   PetscSection   section;

1057:   DMGetLocalSection(dm,&section);
1058:   PetscSectionGetFieldComponents(section,fieldNum,&numComponents);
1059:   DMGetNumFields(dm,&dsNumFields);
1060:   if (dsNumFields) {DMGetField(dm,fieldNum,NULL,&disc);}
1061:   if (disc) {
1062:     PetscObjectGetClassId(disc,&id);
1063:     isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1064:   }
1065:   if (!disc || isContainer) {
1066:     MPI_Comm        comm = PetscObjectComm((PetscObject) dm);
1067:     PetscInt        cStart, cEnd, dim, cellHeight;
1068:     PetscInt        localConeSize = 0, coneSize;
1069:     PetscFE         fe;
1070:     PetscDualSpace  Q;
1071:     PetscSpace      P;
1072:     DM              K;
1073:     PetscQuadrature quad, fquad;
1074:     PetscBool       isSimplex;

1076:     DMPlexGetVTKCellHeight(dm, &cellHeight);
1077:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1078:     DMGetDimension(dm, &dim);
1079:     if (cEnd > cStart) {
1080:       DMPlexGetConeSize(dm, cStart, &localConeSize);
1081:     }
1082:     MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);
1083:     isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE;
1084:     PetscSpaceCreate(comm, &P);
1085:     PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);
1086:     PetscSpaceSetDegree(P, 1, PETSC_DETERMINE);
1087:     PetscSpaceSetNumComponents(P, numComponents);
1088:     PetscSpaceSetNumVariables(P, dim);
1089:     PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);
1090:     PetscSpaceSetUp(P);
1091:     PetscDualSpaceCreate(comm, &Q);
1092:     PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1093:     PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1094:     PetscDualSpaceSetDM(Q, K);
1095:     DMDestroy(&K);
1096:     PetscDualSpaceSetNumComponents(Q, numComponents);
1097:     PetscDualSpaceSetOrder(Q, 1);
1098:     PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);
1099:     PetscDualSpaceSetUp(Q);
1100:     PetscFECreate(comm, &fe);
1101:     PetscFESetType(fe,PETSCFEBASIC);
1102:     PetscFESetBasisSpace(fe, P);
1103:     PetscFESetDualSpace(fe, Q);
1104:     PetscFESetNumComponents(fe, numComponents);
1105:     PetscFESetUp(fe);
1106:     PetscSpaceDestroy(&P);
1107:     PetscDualSpaceDestroy(&Q);
1108:     if (isSimplex) {
1109:       PetscDTGaussJacobiQuadrature(dim,   1, 1, -1.0, 1.0, &quad);
1110:       PetscDTGaussJacobiQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);
1111:     }
1112:     else {
1113:       PetscDTGaussTensorQuadrature(dim,   1, 1, -1.0, 1.0, &quad);
1114:       PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);
1115:     }
1116:     PetscFESetQuadrature(fe, quad);
1117:     PetscFESetFaceQuadrature(fe, fquad);
1118:     PetscQuadratureDestroy(&quad);
1119:     PetscQuadratureDestroy(&fquad);
1120:     disc = (PetscObject) fe;
1121:   } else {
1122:     PetscObjectReference(disc);
1123:   }
1124:   PetscObjectGetClassId(disc,&id);
1125:   if (id == PETSCFE_CLASSID) {
1126:     PetscFE fe = (PetscFE) disc;

1128:     PetscFEGetNumComponents(fe,&numComponents);
1129:   } else {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");}
1130:   DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);
1131:   DMFieldSetType(b,DMFIELDDS);
1132:   dsfield = (DMField_DS *) b->data;
1133:   dsfield->fieldNum = fieldNum;
1134:   DMGetDimension(dm,&dsfield->height);
1135:   dsfield->height++;
1136:   PetscCalloc1(dsfield->height,&dsfield->disc);
1137:   dsfield->disc[0] = disc;
1138:   PetscObjectReference((PetscObject)vec);
1139:   dsfield->vec = vec;
1140:   *field = b;

1142:   return(0);
1143: }