Actual source code: dmfieldds.c

petsc-3.13.6 2020-09-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: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
 96: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
 97: {
 98:   DMField_DS      *dsfield = (DMField_DS *) field->data;
 99:   DM              dm;
100:   PetscObject     disc;
101:   PetscClassId    classid;
102:   PetscInt        nq, nc, dim, meshDim, numCells;
103:   PetscSection    section;
104:   const PetscReal *qpoints;
105:   PetscBool       isStride;
106:   const PetscInt  *points = NULL;
107:   PetscInt        sfirst = -1, stride = -1;
108:   PetscErrorCode  ierr;

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

133:     PetscFEGetDimension(fe,&feDim);
134:     PetscFECreateTabulation(fe,1,nq,qpoints,K,&T);
135:     for (i = 0; i < numCells; i++) {
136:       PetscInt     c = isStride ? (sfirst + i * stride) : points[i];
137:       PetscInt     closureSize;
138:       PetscScalar *elem = NULL;

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

145:           DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,(PetscScalar));
146:         } else {
147:           PetscReal *cB = &((PetscReal *) B)[nc * nq * i];

149:           DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,PetscRealPart);
150:         }
151:       }
152:       if (D) {
153:         if (type == PETSC_SCALAR) {
154:           PetscScalar *cD = &((PetscScalar *) D)[nc * nq * dim * i];

156:           DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),(PetscScalar));
157:         } else {
158:           PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i];

160:           DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),PetscRealPart);
161:         }
162:       }
163:       if (H) {
164:         if (type == PETSC_SCALAR) {
165:           PetscScalar *cH = &((PetscScalar *) H)[nc * nq * dim * dim * i];

167:           DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),(PetscScalar));
168:         } else {
169:           PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i];

171:           DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),PetscRealPart);
172:         }
173:       }
174:       DMPlexVecRestoreClosure(dm,section,dsfield->vec,c,&closureSize,&elem);
175:     }
176:     PetscTabulationDestroy(&T);
177:   } else {SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");}
178:   if (!isStride) {
179:     ISRestoreIndices(pointIS,&points);
180:   }
181:   return(0);
182: }

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

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

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

244:     if (nq) {
245:       PetscInt          K = H ? 2 : (D ? 1 : (B ? 0 : -1));
246:       PetscTabulation   T;
247:       PetscInt     closureSize;
248:       PetscScalar *elem = NULL;
249:       PetscReal   *quadPoints;
250:       PetscQuadrature quad;
251:       PetscInt d, e, f, g;

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

267:           DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,(PetscScalar));
268:         } else {
269:           PetscReal *cB = &cellBr[nc * offset];

271:           DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,PetscRealPart);
272:         }
273:       }
274:       if (D) {
275:         if (datatype == PETSC_SCALAR) {
276:           PetscScalar *cD = &cellDs[nc * dim * offset];

278:           DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),(PetscScalar));
279:           for (p = 0; p < nq; p++) {
280:             for (g = 0; g < nc; g++) {
281:               PetscScalar vs[3];

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

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

317:           DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),(PetscScalar));
318:           for (p = 0; p < nq; p++) {
319:             for (g = 0; g < nc * dimR; g++) {
320:               PetscScalar vs[3];

322:               for (d = 0; d < dimR; d++) {
323:                 vs[d] = 0.;
324:                 for (e = 0; e < dimR; e++) {
325:                   vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
326:                 }
327:               }
328:               for (d = 0; d < dimR; d++) {
329:                 cH[(nc * dimR * p + g) * dimR + d] = vs[d];
330:               }
331:             }
332:             for (g = 0; g < nc; g++) {
333:               for (f = 0; f < dimR; f++) {
334:                 PetscScalar vs[3];

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

351:           DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),PetscRealPart);
352:           for (p = 0; p < nq; p++) {
353:             for (g = 0; g < nc * dimR; g++) {
354:               for (d = 0; d < dimR; d++) {
355:                 v[d] = 0.;
356:                 for (e = 0; e < dimR; e++) {
357:                   v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
358:                 }
359:               }
360:               for (d = 0; d < dimR; d++) {
361:                 cH[(nc * dimR * p + g) * dimR + d] = v[d];
362:               }
363:             }
364:             for (g = 0; g < nc; g++) {
365:               for (f = 0; f < dimR; f++) {
366:                 for (d = 0; d < dimR; d++) {
367:                   v[d] = 0.;
368:                   for (e = 0; e < dimR; e++) {
369:                     v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
370:                   }
371:                 }
372:                 for (d = 0; d < dimR; d++) {
373:                   cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
374:                 }
375:               }
376:             }
377:           }
378:         }
379:       }
380:       DMPlexVecRestoreClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);
381:       PetscTabulationDestroy(&T);
382:     }
383:     offset += nq;
384:   }
385:   {
386:     MPI_Datatype origtype;

388:     if (datatype == PETSC_SCALAR) {
389:       origtype = MPIU_SCALAR;
390:     } else {
391:       origtype = MPIU_REAL;
392:     }
393:     if (B) {
394:       MPI_Datatype Btype;

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

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

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

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

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

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

482:     if (type == PETSC_SCALAR) {
483:       PetscScalar * sB  = (PetscScalar *) B;
484:       PetscScalar * sqB = (PetscScalar *) qB;

486:       for (i = 0; i < numPoints; i++) {
487:         PetscReal vol = 0.;

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

502:       for (i = 0; i < numPoints; i++) {
503:         PetscReal vol = 0.;

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

519:     if (type == PETSC_SCALAR) {
520:       PetscScalar * sD  = (PetscScalar *) D;
521:       PetscScalar * sqD = (PetscScalar *) qD;

523:       for (i = 0; i < numPoints; i++) {
524:         PetscReal vol = 0.;

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

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

548:       for (i = 0; i < numPoints; i++) {
549:         PetscReal vol = 0.;

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

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

574:     if (type == PETSC_SCALAR) {
575:       PetscScalar * sH  = (PetscScalar *) H;
576:       PetscScalar * sqH = (PetscScalar *) qH;

578:       for (i = 0; i < numPoints; i++) {
579:         PetscReal vol = 0.;

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

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

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

612:       for (i = 0; i < numPoints; i++) {
613:         PetscReal vol = 0.;

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

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

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

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

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

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

679:     PetscFEGetBasisSpace(fe, &sp);
680:     PetscSpaceGetDegree(sp, minDegree, maxDegree);
681:   }
682:   return(0);
683: }

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1141:   return(0);
1142: }