Actual source code: dualspace.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petscdmplex.h>

  4: PetscClassId PETSCDUALSPACE_CLASSID = 0;

  6: PetscFunctionList PetscDualSpaceList              = NULL;
  7: PetscBool         PetscDualSpaceRegisterAllCalled = PETSC_FALSE;

  9: /*@C
 10:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

 12:   Not Collective

 14:   Input Parameters:
 15: + name        - The name of a new user-defined creation routine
 16: - create_func - The creation routine itself

 18:   Notes:
 19:   PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces

 21:   Sample usage:
 22: .vb
 23:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
 24: .ve

 26:   Then, your PetscDualSpace type can be chosen with the procedural interface via
 27: .vb
 28:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
 29:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
 30: .ve
 31:    or at runtime via the option
 32: .vb
 33:     -petscdualspace_type my_dual_space
 34: .ve

 36:   Level: advanced

 38: .keywords: PetscDualSpace, register
 39: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()

 41: @*/
 42: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
 43: {

 47:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
 48:   return(0);
 49: }

 51: /*@C
 52:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

 54:   Collective on PetscDualSpace

 56:   Input Parameters:
 57: + sp   - The PetscDualSpace object
 58: - name - The kind of space

 60:   Options Database Key:
 61: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types

 63:   Level: intermediate

 65: .keywords: PetscDualSpace, set, type
 66: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
 67: @*/
 68: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
 69: {
 70:   PetscErrorCode (*r)(PetscDualSpace);
 71:   PetscBool      match;

 76:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
 77:   if (match) return(0);

 79:   if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
 80:   PetscFunctionListFind(PetscDualSpaceList, name, &r);
 81:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);

 83:   if (sp->ops->destroy) {
 84:     (*sp->ops->destroy)(sp);
 85:     sp->ops->destroy = NULL;
 86:   }
 87:   (*r)(sp);
 88:   PetscObjectChangeTypeName((PetscObject) sp, name);
 89:   return(0);
 90: }

 92: /*@C
 93:   PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.

 95:   Not Collective

 97:   Input Parameter:
 98: . sp  - The PetscDualSpace

100:   Output Parameter:
101: . name - The PetscDualSpace type name

103:   Level: intermediate

105: .keywords: PetscDualSpace, get, type, name
106: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
107: @*/
108: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
109: {

115:   if (!PetscDualSpaceRegisterAllCalled) {
116:     PetscDualSpaceRegisterAll();
117:   }
118:   *name = ((PetscObject) sp)->type_name;
119:   return(0);
120: }

122: /*@
123:   PetscDualSpaceView - Views a PetscDualSpace

125:   Collective on PetscDualSpace

127:   Input Parameter:
128: + sp - the PetscDualSpace object to view
129: - v  - the viewer

131:   Level: developer

133: .seealso PetscDualSpaceDestroy()
134: @*/
135: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
136: {
137:   PetscBool      iascii;

143:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
144:   PetscObjectPrintClassNamePrefixType((PetscObject)sp, v);
145:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
146:   PetscViewerASCIIPushTab(v);
147:   if (iascii) {PetscViewerASCIIPrintf(v, "Dual space of order %D with %D components\n", sp->order, sp->Nc);}
148:   if (sp->ops->view) {(*sp->ops->view)(sp, v);}
149:   PetscViewerASCIIPopTab(v);
150:   return(0);
151: }

153: /*@
154:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

156:   Collective on PetscDualSpace

158:   Input Parameter:
159: . sp - the PetscDualSpace object to set options for

161:   Options Database:
162: . -petscspace_degree the approximation order of the space

164:   Level: developer

166: .seealso PetscDualSpaceView()
167: @*/
168: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
169: {
170:   const char    *defaultType;
171:   char           name[256];
172:   PetscBool      flg;

177:   if (!((PetscObject) sp)->type_name) {
178:     defaultType = PETSCDUALSPACELAGRANGE;
179:   } else {
180:     defaultType = ((PetscObject) sp)->type_name;
181:   }
182:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

184:   PetscObjectOptionsBegin((PetscObject) sp);
185:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
186:   if (flg) {
187:     PetscDualSpaceSetType(sp, name);
188:   } else if (!((PetscObject) sp)->type_name) {
189:     PetscDualSpaceSetType(sp, defaultType);
190:   }
191:   PetscOptionsInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
192:   PetscOptionsInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL);
193:   if (sp->ops->setfromoptions) {
194:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
195:   }
196:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
197:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
198:   PetscOptionsEnd();
199:   PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
200:   return(0);
201: }

203: /*@
204:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

206:   Collective on PetscDualSpace

208:   Input Parameter:
209: . sp - the PetscDualSpace object to setup

211:   Level: developer

213: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
214: @*/
215: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
216: {

221:   if (sp->setupcalled) return(0);
222:   sp->setupcalled = PETSC_TRUE;
223:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
224:   return(0);
225: }

227: /*@
228:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

230:   Collective on PetscDualSpace

232:   Input Parameter:
233: . sp - the PetscDualSpace object to destroy

235:   Level: developer

237: .seealso PetscDualSpaceView()
238: @*/
239: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
240: {
241:   PetscInt       dim, f;

245:   if (!*sp) return(0);

248:   if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
249:   ((PetscObject) (*sp))->refct = 0;

251:   PetscDualSpaceGetDimension(*sp, &dim);
252:   for (f = 0; f < dim; ++f) {
253:     PetscQuadratureDestroy(&(*sp)->functional[f]);
254:   }
255:   PetscFree((*sp)->functional);
256:   PetscQuadratureDestroy(&(*sp)->allPoints);
257:   DMDestroy(&(*sp)->dm);

259:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
260:   PetscHeaderDestroy(sp);
261:   return(0);
262: }

264: /*@
265:   PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().

267:   Collective on MPI_Comm

269:   Input Parameter:
270: . comm - The communicator for the PetscDualSpace object

272:   Output Parameter:
273: . sp - The PetscDualSpace object

275:   Level: beginner

277: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
278: @*/
279: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
280: {
281:   PetscDualSpace s;

286:   PetscCitationsRegister(FECitation,&FEcite);
287:   *sp  = NULL;
288:   PetscFEInitializePackage();

290:   PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);

292:   s->order = 0;
293:   s->Nc    = 1;
294:   s->setupcalled = PETSC_FALSE;

296:   *sp = s;
297:   return(0);
298: }

300: /*@
301:   PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.

303:   Collective on PetscDualSpace

305:   Input Parameter:
306: . sp - The original PetscDualSpace

308:   Output Parameter:
309: . spNew - The duplicate PetscDualSpace

311:   Level: beginner

313: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
314: @*/
315: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
316: {

322:   (*sp->ops->duplicate)(sp, spNew);
323:   return(0);
324: }

326: /*@
327:   PetscDualSpaceGetDM - Get the DM representing the reference cell

329:   Not collective

331:   Input Parameter:
332: . sp - The PetscDualSpace

334:   Output Parameter:
335: . dm - The reference cell

337:   Level: intermediate

339: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
340: @*/
341: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
342: {
346:   *dm = sp->dm;
347:   return(0);
348: }

350: /*@
351:   PetscDualSpaceSetDM - Get the DM representing the reference cell

353:   Not collective

355:   Input Parameters:
356: + sp - The PetscDualSpace
357: - dm - The reference cell

359:   Level: intermediate

361: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
362: @*/
363: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
364: {

370:   DMDestroy(&sp->dm);
371:   PetscObjectReference((PetscObject) dm);
372:   sp->dm = dm;
373:   return(0);
374: }

376: /*@
377:   PetscDualSpaceGetOrder - Get the order of the dual space

379:   Not collective

381:   Input Parameter:
382: . sp - The PetscDualSpace

384:   Output Parameter:
385: . order - The order

387:   Level: intermediate

389: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
390: @*/
391: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
392: {
396:   *order = sp->order;
397:   return(0);
398: }

400: /*@
401:   PetscDualSpaceSetOrder - Set the order of the dual space

403:   Not collective

405:   Input Parameters:
406: + sp - The PetscDualSpace
407: - order - The order

409:   Level: intermediate

411: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
412: @*/
413: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
414: {
417:   sp->order = order;
418:   return(0);
419: }

421: /*@
422:   PetscDualSpaceGetNumComponents - Return the number of components for this space

424:   Input Parameter:
425: . sp - The PetscDualSpace

427:   Output Parameter:
428: . Nc - The number of components

430:   Note: A vector space, for example, will have d components, where d is the spatial dimension

432:   Level: intermediate

434: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
435: @*/
436: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
437: {
441:   *Nc = sp->Nc;
442:   return(0);
443: }

445: /*@
446:   PetscDualSpaceSetNumComponents - Set the number of components for this space

448:   Input Parameters:
449: + sp - The PetscDualSpace
450: - order - The number of components

452:   Level: intermediate

454: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
455: @*/
456: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
457: {
460:   sp->Nc = Nc;
461:   return(0);
462: }

464: /*@
465:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

467:   Not collective

469:   Input Parameter:
470: . sp - The PetscDualSpace

472:   Output Parameter:
473: . tensor - Whether the dual space has tensor layout (vs. simplicial)

475:   Level: intermediate

477: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
478: @*/
479: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
480: {

486:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
487:   return(0);
488: }

490: /*@
491:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

493:   Not collective

495:   Input Parameters:
496: + sp - The PetscDualSpace
497: - tensor - Whether the dual space has tensor layout (vs. simplicial)

499:   Level: intermediate

501: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
502: @*/
503: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
504: {

509:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
510:   return(0);
511: }

513: /*@
514:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

516:   Not collective

518:   Input Parameters:
519: + sp - The PetscDualSpace
520: - i  - The basis number

522:   Output Parameter:
523: . functional - The basis functional

525:   Level: intermediate

527: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
528: @*/
529: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
530: {
531:   PetscInt       dim;

537:   PetscDualSpaceGetDimension(sp, &dim);
538:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
539:   *functional = sp->functional[i];
540:   return(0);
541: }

543: /*@
544:   PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals

546:   Not collective

548:   Input Parameter:
549: . sp - The PetscDualSpace

551:   Output Parameter:
552: . dim - The dimension

554:   Level: intermediate

556: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
557: @*/
558: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
559: {

565:   *dim = 0;
566:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
567:   return(0);
568: }

570: /*@C
571:   PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension

573:   Not collective

575:   Input Parameter:
576: . sp - The PetscDualSpace

578:   Output Parameter:
579: . numDof - An array of length dim+1 which holds the number of dofs for each dimension

581:   Level: intermediate

583: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
584: @*/
585: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
586: {

592:   (*sp->ops->getnumdof)(sp, numDof);
593:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
594:   return(0);
595: }

597: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
598: {
599:   DM             dm;
600:   PetscInt       pStart, pEnd, depth, h, offset;
601:   const PetscInt *numDof;

605:   PetscDualSpaceGetDM(sp,&dm);
606:   DMPlexGetChart(dm,&pStart,&pEnd);
607:   PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
608:   PetscSectionSetChart(*section,pStart,pEnd);
609:   DMPlexGetDepth(dm,&depth);
610:   PetscDualSpaceGetNumDof(sp,&numDof);
611:   for (h = 0; h <= depth; h++) {
612:     PetscInt hStart, hEnd, p, dof;

614:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
615:     dof = numDof[depth - h];
616:     for (p = hStart; p < hEnd; p++) {
617:       PetscSectionSetDof(*section,p,dof);
618:     }
619:   }
620:   PetscSectionSetUp(*section);
621:   for (h = 0, offset = 0; h <= depth; h++) {
622:     PetscInt hStart, hEnd, p, dof;

624:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
625:     dof = numDof[depth - h];
626:     for (p = hStart; p < hEnd; p++) {
627:       PetscSectionGetDof(*section,p,&dof);
628:       PetscSectionSetOffset(*section,p,offset);
629:       offset += dof;
630:     }
631:   }
632:   return(0);
633: }

635: /*@
636:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

638:   Collective on PetscDualSpace

640:   Input Parameters:
641: + sp      - The PetscDualSpace
642: . dim     - The spatial dimension
643: - simplex - Flag for simplex, otherwise use a tensor-product cell

645:   Output Parameter:
646: . refdm - The reference cell

648:   Level: advanced

650: .keywords: PetscDualSpace, reference cell
651: .seealso: PetscDualSpaceCreate(), DMPLEX
652: @*/
653: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
654: {

658:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
659:   return(0);
660: }

662: /*@C
663:   PetscDualSpaceApply - Apply a functional from the dual space basis to an input function

665:   Input Parameters:
666: + sp      - The PetscDualSpace object
667: . f       - The basis functional index
668: . time    - The time
669: . cgeom   - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian) (or evaluated at the coordinates of the functional)
670: . numComp - The number of components for the function
671: . func    - The input function
672: - ctx     - A context for the function

674:   Output Parameter:
675: . value   - numComp output values

677:   Note: The calling sequence for the callback func is given by:

679: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
680: $      PetscInt numComponents, PetscScalar values[], void *ctx)

682:   Level: developer

684: .seealso: PetscDualSpaceCreate()
685: @*/
686: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
687: {

694:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
695:   return(0);
696: }

698: /*@C
699:   PetscDualSpaceApplyAll - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()

701:   Input Parameters:
702: + sp        - The PetscDualSpace object
703: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

705:   Output Parameter:
706: . spValue   - The values of all dual space functionals

708:   Level: developer

710: .seealso: PetscDualSpaceCreate()
711: @*/
712: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
713: {

718:   (*sp->ops->applyall)(sp, pointEval, spValue);
719:   return(0);
720: }

722: /*@C
723:   PetscDualSpaceApplyDefault - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional.

725:   Input Parameters:
726: + sp    - The PetscDualSpace object
727: . f     - The basis functional index
728: . time  - The time
729: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
730: . Nc    - The number of components for the function
731: . func  - The input function
732: - ctx   - A context for the function

734:   Output Parameter:
735: . value   - The output value

737:   Note: The calling sequence for the callback func is given by:

739: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
740: $      PetscInt numComponents, PetscScalar values[], void *ctx)

742: and the idea is to evaluate the functional as an integral

744: $ n(f) = int dx n(x) . f(x)

746: where both n and f have Nc components.

748:   Level: developer

750: .seealso: PetscDualSpaceCreate()
751: @*/
752: PetscErrorCode PetscDualSpaceApplyDefault(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFEGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
753: {
754:   DM               dm;
755:   PetscQuadrature  n;
756:   const PetscReal *points, *weights;
757:   PetscReal        x[3];
758:   PetscScalar     *val;
759:   PetscInt         dim, dE, qNc, c, Nq, q;
760:   PetscBool        isAffine;
761:   PetscErrorCode   ierr;

766:   PetscDualSpaceGetDM(sp, &dm);
767:   PetscDualSpaceGetFunctional(sp, f, &n);
768:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
769:   if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
770:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
771:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
772:   *value = 0.0;
773:   isAffine = cgeom->isAffine;
774:   dE = cgeom->dimEmbed;
775:   for (q = 0; q < Nq; ++q) {
776:     if (isAffine) {
777:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
778:       (*func)(dE, time, x, Nc, val, ctx);
779:     } else {
780:       (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
781:     }
782:     for (c = 0; c < Nc; ++c) {
783:       *value += val[c]*weights[q*Nc+c];
784:     }
785:   }
786:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
787:   return(0);
788: }

790: /*@C
791:   PetscDualSpaceApplyAllDefault - Apply all functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetAllPoints()

793:   Input Parameters:
794: + sp        - The PetscDualSpace object
795: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

797:   Output Parameter:
798: . spValue   - The values of all dual space functionals

800:   Level: developer

802: .seealso: PetscDualSpaceCreate()
803: @*/
804: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
805: {
806:   PetscQuadrature  n;
807:   const PetscReal *points, *weights;
808:   PetscInt         qNc, c, Nq, q, f, spdim, Nc;
809:   PetscInt         offset;
810:   PetscErrorCode   ierr;

816:   PetscDualSpaceGetDimension(sp, &spdim);
817:   PetscDualSpaceGetNumComponents(sp, &Nc);
818:   for (f = 0, offset = 0; f < spdim; f++) {
819:     PetscDualSpaceGetFunctional(sp, f, &n);
820:     PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
821:     if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
822:     spValue[f] = 0.0;
823:     for (q = 0; q < Nq; ++q) {
824:       for (c = 0; c < Nc; ++c) {
825:         spValue[f] += pointEval[offset++]*weights[q*Nc+c];
826:       }
827:     }
828:   }
829:   return(0);
830: }

832: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
833: {

839:   if (!sp->allPoints && sp->ops->createallpoints) {
840:     (*sp->ops->createallpoints)(sp,&sp->allPoints);
841:   }
842:   *allPoints = sp->allPoints;
843:   return(0);
844: }

846: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
847: {
848:   PetscInt        spdim;
849:   PetscInt        numPoints, offset;
850:   PetscReal       *points;
851:   PetscInt        f, dim;
852:   PetscQuadrature q;
853:   PetscErrorCode  ierr;

856:   PetscDualSpaceGetDimension(sp,&spdim);
857:   if (!spdim) {
858:     PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
859:     PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
860:   }
861:   PetscDualSpaceGetFunctional(sp,0,&q);
862:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
863:   for (f = 1; f < spdim; f++) {
864:     PetscInt Np;

866:     PetscDualSpaceGetFunctional(sp,f,&q);
867:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
868:     numPoints += Np;
869:   }
870:   PetscMalloc1(dim*numPoints,&points);
871:   for (f = 0, offset = 0; f < spdim; f++) {
872:     const PetscReal *p;
873:     PetscInt        Np, i;

875:     PetscDualSpaceGetFunctional(sp,f,&q);
876:     PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
877:     for (i = 0; i < Np * dim; i++) {
878:       points[offset + i] = p[i];
879:     }
880:     offset += Np * dim;
881:   }
882:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
883:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
884:   return(0);
885: }

887: /*@C
888:   PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function by assuming a point evaluation functional at the cell centroid.

890:   Input Parameters:
891: + sp    - The PetscDualSpace object
892: . f     - The basis functional index
893: . time  - The time
894: . cgeom - A context with geometric information for this cell, we currently just use the centroid
895: . Nc    - The number of components for the function
896: . func  - The input function
897: - ctx   - A context for the function

899:   Output Parameter:
900: . value - The output value (scalar)

902:   Note: The calling sequence for the callback func is given by:

904: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
905: $      PetscInt numComponents, PetscScalar values[], void *ctx)

907: and the idea is to evaluate the functional as an integral

909: $ n(f) = int dx n(x) . f(x)

911: where both n and f have Nc components.

913:   Level: developer

915: .seealso: PetscDualSpaceCreate()
916: @*/
917: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *cgeom, PetscInt Nc, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
918: {
919:   DM               dm;
920:   PetscQuadrature  n;
921:   const PetscReal *points, *weights;
922:   PetscScalar     *val;
923:   PetscInt         dimEmbed, qNc, c, Nq, q;
924:   PetscErrorCode   ierr;

929:   PetscDualSpaceGetDM(sp, &dm);
930:   DMGetCoordinateDim(dm, &dimEmbed);
931:   PetscDualSpaceGetFunctional(sp, f, &n);
932:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
933:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
934:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
935:   *value = 0.;
936:   for (q = 0; q < Nq; ++q) {
937:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
938:     for (c = 0; c < Nc; ++c) {
939:       *value += val[c]*weights[q*Nc+c];
940:     }
941:   }
942:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
943:   return(0);
944: }

946: /*@
947:   PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a
948:   given height.  This assumes that the reference cell is symmetric over points of this height.

950:   If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
951:   pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
952:   support extracting subspaces, then NULL is returned.

954:   This does not increment the reference count on the returned dual space, and the user should not destroy it.

956:   Not collective

958:   Input Parameters:
959: + sp - the PetscDualSpace object
960: - height - the height of the mesh point for which the subspace is desired

962:   Output Parameter:
963: . subsp - the subspace.  Note that the functionals in the subspace are with respect to the intrinsic geometry of the
964:   point, which will be of lesser dimension if height > 0.

966:   Level: advanced

968: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
969: @*/
970: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
971: {

977:   *subsp = NULL;
978:   if (sp->ops->getheightsubspace) {
979:     (*sp->ops->getheightsubspace)(sp, height, subsp);
980:   }
981:   return(0);
982: }

984: /*@
985:   PetscDualSpaceGetPointSubspace - Get the subset of the dual space basis that is supported on a particular mesh point.

987:   If the dual space is not defined on the mesh point (e.g. if the space is discontinuous and pointwise values are not
988:   defined on the element boundaries), or if the implementation of PetscDualSpace does not support extracting
989:   subspaces, then NULL is returned.

991:   This does not increment the reference count on the returned dual space, and the user should not destroy it.

993:   Not collective

995:   Input Parameters:
996: + sp - the PetscDualSpace object
997: - point - the point (in the dual space's DM) for which the subspace is desired

999:   Output Parameters:
1000:   bdsp - the subspace.  Note that the functionals in the subspace are with respect to the intrinsic geometry of the
1001:   point, which will be of lesser dimension if height > 0.

1003:   Level: advanced

1005: .seealso: PetscDualSpace
1006: @*/
1007: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1008: {

1014:   *bdsp = NULL;
1015:   if (sp->ops->getpointsubspace) {
1016:     (*sp->ops->getpointsubspace)(sp,point,bdsp);
1017:   } else if (sp->ops->getheightsubspace) {
1018:     DM       dm;
1019:     DMLabel  label;
1020:     PetscInt dim, depth, height;

1022:     PetscDualSpaceGetDM(sp,&dm);
1023:     DMPlexGetDepth(dm,&dim);
1024:     DMPlexGetDepthLabel(dm,&label);
1025:     DMLabelGetValue(label,point,&depth);
1026:     height = dim - depth;
1027:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1028:   }
1029:   return(0);
1030: }