Actual source code: dualspace.c

petsc-3.10.2 2018-10-09
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: {

141:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
142:   if (sp->ops->view) {(*sp->ops->view)(sp, v);}
143:   return(0);
144: }

146: /*@
147:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

149:   Collective on PetscDualSpace

151:   Input Parameter:
152: . sp - the PetscDualSpace object to set options for

154:   Options Database:
155: . -petscspace_degree the approximation order of the space

157:   Level: developer

159: .seealso PetscDualSpaceView()
160: @*/
161: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
162: {
163:   const char    *defaultType;
164:   char           name[256];
165:   PetscBool      flg;

170:   if (!((PetscObject) sp)->type_name) {
171:     defaultType = PETSCDUALSPACELAGRANGE;
172:   } else {
173:     defaultType = ((PetscObject) sp)->type_name;
174:   }
175:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

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

196: /*@
197:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

199:   Collective on PetscDualSpace

201:   Input Parameter:
202: . sp - the PetscDualSpace object to setup

204:   Level: developer

206: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
207: @*/
208: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
209: {

214:   if (sp->setupcalled) return(0);
215:   sp->setupcalled = PETSC_TRUE;
216:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
217:   return(0);
218: }

220: /*@
221:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

223:   Collective on PetscDualSpace

225:   Input Parameter:
226: . sp - the PetscDualSpace object to destroy

228:   Level: developer

230: .seealso PetscDualSpaceView()
231: @*/
232: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
233: {
234:   PetscInt       dim, f;

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

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

244:   PetscDualSpaceGetDimension(*sp, &dim);
245:   for (f = 0; f < dim; ++f) {
246:     PetscQuadratureDestroy(&(*sp)->functional[f]);
247:   }
248:   PetscFree((*sp)->functional);
249:   PetscQuadratureDestroy(&(*sp)->allPoints);
250:   DMDestroy(&(*sp)->dm);

252:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
253:   PetscHeaderDestroy(sp);
254:   return(0);
255: }

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

260:   Collective on MPI_Comm

262:   Input Parameter:
263: . comm - The communicator for the PetscDualSpace object

265:   Output Parameter:
266: . sp - The PetscDualSpace object

268:   Level: beginner

270: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
271: @*/
272: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
273: {
274:   PetscDualSpace s;

279:   PetscCitationsRegister(FECitation,&FEcite);
280:   *sp  = NULL;
281:   PetscFEInitializePackage();

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

285:   s->order = 0;
286:   s->Nc    = 1;
287:   s->setupcalled = PETSC_FALSE;

289:   *sp = s;
290:   return(0);
291: }

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

296:   Collective on PetscDualSpace

298:   Input Parameter:
299: . sp - The original PetscDualSpace

301:   Output Parameter:
302: . spNew - The duplicate PetscDualSpace

304:   Level: beginner

306: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
307: @*/
308: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
309: {

315:   (*sp->ops->duplicate)(sp, spNew);
316:   return(0);
317: }

319: /*@
320:   PetscDualSpaceGetDM - Get the DM representing the reference cell

322:   Not collective

324:   Input Parameter:
325: . sp - The PetscDualSpace

327:   Output Parameter:
328: . dm - The reference cell

330:   Level: intermediate

332: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
333: @*/
334: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
335: {
339:   *dm = sp->dm;
340:   return(0);
341: }

343: /*@
344:   PetscDualSpaceSetDM - Get the DM representing the reference cell

346:   Not collective

348:   Input Parameters:
349: + sp - The PetscDualSpace
350: - dm - The reference cell

352:   Level: intermediate

354: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
355: @*/
356: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
357: {

363:   DMDestroy(&sp->dm);
364:   PetscObjectReference((PetscObject) dm);
365:   sp->dm = dm;
366:   return(0);
367: }

369: /*@
370:   PetscDualSpaceGetOrder - Get the order of the dual space

372:   Not collective

374:   Input Parameter:
375: . sp - The PetscDualSpace

377:   Output Parameter:
378: . order - The order

380:   Level: intermediate

382: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
383: @*/
384: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
385: {
389:   *order = sp->order;
390:   return(0);
391: }

393: /*@
394:   PetscDualSpaceSetOrder - Set the order of the dual space

396:   Not collective

398:   Input Parameters:
399: + sp - The PetscDualSpace
400: - order - The order

402:   Level: intermediate

404: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
405: @*/
406: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
407: {
410:   sp->order = order;
411:   return(0);
412: }

414: /*@
415:   PetscDualSpaceGetNumComponents - Return the number of components for this space

417:   Input Parameter:
418: . sp - The PetscDualSpace

420:   Output Parameter:
421: . Nc - The number of components

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

425:   Level: intermediate

427: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
428: @*/
429: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
430: {
434:   *Nc = sp->Nc;
435:   return(0);
436: }

438: /*@
439:   PetscDualSpaceSetNumComponents - Set the number of components for this space

441:   Input Parameters:
442: + sp - The PetscDualSpace
443: - order - The number of components

445:   Level: intermediate

447: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
448: @*/
449: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
450: {
453:   sp->Nc = Nc;
454:   return(0);
455: }

457: /*@
458:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

460:   Not collective

462:   Input Parameter:
463: . sp - The PetscDualSpace

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

468:   Level: intermediate

470: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
471: @*/
472: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
473: {

479:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
480:   return(0);
481: }

483: /*@
484:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

486:   Not collective

488:   Input Parameters:
489: + sp - The PetscDualSpace
490: - tensor - Whether the dual space has tensor layout (vs. simplicial)

492:   Level: intermediate

494: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
495: @*/
496: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
497: {

502:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
503:   return(0);
504: }

506: /*@
507:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

509:   Not collective

511:   Input Parameters:
512: + sp - The PetscDualSpace
513: - i  - The basis number

515:   Output Parameter:
516: . functional - The basis functional

518:   Level: intermediate

520: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
521: @*/
522: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
523: {
524:   PetscInt       dim;

530:   PetscDualSpaceGetDimension(sp, &dim);
531:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
532:   *functional = sp->functional[i];
533:   return(0);
534: }

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

539:   Not collective

541:   Input Parameter:
542: . sp - The PetscDualSpace

544:   Output Parameter:
545: . dim - The dimension

547:   Level: intermediate

549: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
550: @*/
551: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
552: {

558:   *dim = 0;
559:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
560:   return(0);
561: }

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

566:   Not collective

568:   Input Parameter:
569: . sp - The PetscDualSpace

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

574:   Level: intermediate

576: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
577: @*/
578: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
579: {

585:   (*sp->ops->getnumdof)(sp, numDof);
586:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
587:   return(0);
588: }

590: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
591: {
592:   DM             dm;
593:   PetscInt       pStart, pEnd, depth, h, offset;
594:   const PetscInt *numDof;

598:   PetscDualSpaceGetDM(sp,&dm);
599:   DMPlexGetChart(dm,&pStart,&pEnd);
600:   PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
601:   PetscSectionSetChart(*section,pStart,pEnd);
602:   DMPlexGetDepth(dm,&depth);
603:   PetscDualSpaceGetNumDof(sp,&numDof);
604:   for (h = 0; h <= depth; h++) {
605:     PetscInt hStart, hEnd, p, dof;

607:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
608:     dof = numDof[depth - h];
609:     for (p = hStart; p < hEnd; p++) {
610:       PetscSectionSetDof(*section,p,dof);
611:     }
612:   }
613:   PetscSectionSetUp(*section);
614:   for (h = 0, offset = 0; h <= depth; h++) {
615:     PetscInt hStart, hEnd, p, dof;

617:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
618:     dof = numDof[depth - h];
619:     for (p = hStart; p < hEnd; p++) {
620:       PetscSectionGetDof(*section,p,&dof);
621:       PetscSectionSetOffset(*section,p,offset);
622:       offset += dof;
623:     }
624:   }
625:   return(0);
626: }

628: /*@
629:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

631:   Collective on PetscDualSpace

633:   Input Parameters:
634: + sp      - The PetscDualSpace
635: . dim     - The spatial dimension
636: - simplex - Flag for simplex, otherwise use a tensor-product cell

638:   Output Parameter:
639: . refdm - The reference cell

641:   Level: advanced

643: .keywords: PetscDualSpace, reference cell
644: .seealso: PetscDualSpaceCreate(), DMPLEX
645: @*/
646: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
647: {

651:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
652:   return(0);
653: }

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

658:   Input Parameters:
659: + sp      - The PetscDualSpace object
660: . f       - The basis functional index
661: . time    - The time
662: . 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)
663: . numComp - The number of components for the function
664: . func    - The input function
665: - ctx     - A context for the function

667:   Output Parameter:
668: . value   - numComp output values

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

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

675:   Level: developer

677: .seealso: PetscDualSpaceCreate()
678: @*/
679: 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)
680: {

687:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
688:   return(0);
689: }

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

694:   Input Parameters:
695: + sp        - The PetscDualSpace object
696: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

698:   Output Parameter:
699: . spValue   - The values of all dual space functionals

701:   Level: developer

703: .seealso: PetscDualSpaceCreate()
704: @*/
705: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
706: {

711:   (*sp->ops->applyall)(sp, pointEval, spValue);
712:   return(0);
713: }

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

718:   Input Parameters:
719: + sp    - The PetscDualSpace object
720: . f     - The basis functional index
721: . time  - The time
722: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
723: . Nc    - The number of components for the function
724: . func  - The input function
725: - ctx   - A context for the function

727:   Output Parameter:
728: . value   - The output value

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

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

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

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

739: where both n and f have Nc components.

741:   Level: developer

743: .seealso: PetscDualSpaceCreate()
744: @*/
745: 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)
746: {
747:   DM               dm;
748:   PetscQuadrature  n;
749:   const PetscReal *points, *weights;
750:   PetscReal        x[3];
751:   PetscScalar     *val;
752:   PetscInt         dim, dE, qNc, c, Nq, q;
753:   PetscBool        isAffine;
754:   PetscErrorCode   ierr;

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

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

786:   Input Parameters:
787: + sp        - The PetscDualSpace object
788: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

790:   Output Parameter:
791: . spValue   - The values of all dual space functionals

793:   Level: developer

795: .seealso: PetscDualSpaceCreate()
796: @*/
797: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
798: {
799:   PetscQuadrature  n;
800:   const PetscReal *points, *weights;
801:   PetscInt         qNc, c, Nq, q, f, spdim, Nc;
802:   PetscInt         offset;
803:   PetscErrorCode   ierr;

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

825: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
826: {

832:   if (!sp->allPoints && sp->ops->createallpoints) {
833:     (*sp->ops->createallpoints)(sp,&sp->allPoints);
834:   }
835:   *allPoints = sp->allPoints;
836:   return(0);
837: }

839: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
840: {
841:   PetscInt        spdim;
842:   PetscInt        numPoints, offset;
843:   PetscReal       *points;
844:   PetscInt        f, dim;
845:   PetscQuadrature q;
846:   PetscErrorCode  ierr;

849:   PetscDualSpaceGetDimension(sp,&spdim);
850:   if (!spdim) {
851:     PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
852:     PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
853:   }
854:   PetscDualSpaceGetFunctional(sp,0,&q);
855:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
856:   for (f = 1; f < spdim; f++) {
857:     PetscInt Np;

859:     PetscDualSpaceGetFunctional(sp,f,&q);
860:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
861:     numPoints += Np;
862:   }
863:   PetscMalloc1(dim*numPoints,&points);
864:   for (f = 0, offset = 0; f < spdim; f++) {
865:     const PetscReal *p;
866:     PetscInt        Np, i;

868:     PetscDualSpaceGetFunctional(sp,f,&q);
869:     PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
870:     for (i = 0; i < Np * dim; i++) {
871:       points[offset + i] = p[i];
872:     }
873:     offset += Np * dim;
874:   }
875:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
876:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
877:   return(0);
878: }

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

883:   Input Parameters:
884: + sp    - The PetscDualSpace object
885: . f     - The basis functional index
886: . time  - The time
887: . cgeom - A context with geometric information for this cell, we currently just use the centroid
888: . Nc    - The number of components for the function
889: . func  - The input function
890: - ctx   - A context for the function

892:   Output Parameter:
893: . value - The output value (scalar)

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

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

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

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

904: where both n and f have Nc components.

906:   Level: developer

908: .seealso: PetscDualSpaceCreate()
909: @*/
910: 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)
911: {
912:   DM               dm;
913:   PetscQuadrature  n;
914:   const PetscReal *points, *weights;
915:   PetscScalar     *val;
916:   PetscInt         dimEmbed, qNc, c, Nq, q;
917:   PetscErrorCode   ierr;

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

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

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

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

949:   Not collective

951:   Input Parameters:
952: + sp - the PetscDualSpace object
953: - height - the height of the mesh point for which the subspace is desired

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

959:   Level: advanced

961: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
962: @*/
963: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
964: {

970:   *subsp = NULL;
971:   if (sp->ops->getheightsubspace) {
972:     (*sp->ops->getheightsubspace)(sp, height, subsp);
973:   }
974:   return(0);
975: }

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

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

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

986:   Not collective

988:   Input Parameters:
989: + sp - the PetscDualSpace object
990: - point - the point (in the dual space's DM) for which the subspace is desired

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

996:   Level: advanced

998: .seealso: PetscDualSpace
999: @*/
1000: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1001: {

1007:   *bdsp = NULL;
1008:   if (sp->ops->getpointsubspace) {
1009:     (*sp->ops->getpointsubspace)(sp,point,bdsp);
1010:   } else if (sp->ops->getheightsubspace) {
1011:     DM       dm;
1012:     DMLabel  label;
1013:     PetscInt dim, depth, height;

1015:     PetscDualSpaceGetDM(sp,&dm);
1016:     DMPlexGetDepth(dm,&dim);
1017:     DMPlexGetDepthLabel(dm,&label);
1018:     DMLabelGetValue(label,point,&depth);
1019:     height = dim - depth;
1020:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1021:   }
1022:   return(0);
1023: }