Actual source code: dualspace.c

petsc-master 2019-06-15
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: const char *const PetscDualSpaceReferenceCells[] = {"SIMPLEX", "TENSOR", "PetscDualSpaceReferenceCell", "PETSCDUALSPACE_REFCELL_",0};

 11: /*
 12:   PetscDualSpaceLatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
 13:                                                      Ordering is lexicographic with lowest index as least significant in ordering.
 14:                                                      e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,0}.

 16:   Input Parameters:
 17: + len - The length of the tuple
 18: . max - The maximum sum
 19: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

 21:   Output Parameter:
 22: . tup - A tuple of len integers whos sum is at most 'max'

 24:   Level: developer

 26: .seealso: PetscDualSpaceTensorPointLexicographic_Internal()
 27: */
 28: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
 29: {
 31:   while (len--) {
 32:     max -= tup[len];
 33:     if (!max) {
 34:       tup[len] = 0;
 35:       break;
 36:     }
 37:   }
 38:   tup[++len]++;
 39:   return(0);
 40: }

 42: /*
 43:   PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
 44:                                                     Ordering is lexicographic with lowest index as least significant in ordering.
 45:                                                     e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.

 47:   Input Parameters:
 48: + len - The length of the tuple
 49: . max - The maximum value
 50: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

 52:   Output Parameter:
 53: . tup - A tuple of len integers whos sum is at most 'max'

 55:   Level: developer

 57: .seealso: PetscDualSpaceLatticePointLexicographic_Internal()
 58: */
 59: PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
 60: {
 61:   PetscInt       i;

 64:   for (i = 0; i < len; i++) {
 65:     if (tup[i] < max) {
 66:       break;
 67:     } else {
 68:       tup[i] = 0;
 69:     }
 70:   }
 71:   tup[i]++;
 72:   return(0);
 73: }

 75: /*@C
 76:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

 78:   Not Collective

 80:   Input Parameters:
 81: + name        - The name of a new user-defined creation routine
 82: - create_func - The creation routine itself

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

 87:   Sample usage:
 88: .vb
 89:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
 90: .ve

 92:   Then, your PetscDualSpace type can be chosen with the procedural interface via
 93: .vb
 94:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
 95:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
 96: .ve
 97:    or at runtime via the option
 98: .vb
 99:     -petscdualspace_type my_dual_space
100: .ve

102:   Level: advanced

104: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()

106: @*/
107: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
108: {

112:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
113:   return(0);
114: }

116: /*@C
117:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

119:   Collective on sp

121:   Input Parameters:
122: + sp   - The PetscDualSpace object
123: - name - The kind of space

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

128:   Level: intermediate

130: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
131: @*/
132: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
133: {
134:   PetscErrorCode (*r)(PetscDualSpace);
135:   PetscBool      match;

140:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
141:   if (match) return(0);

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

147:   if (sp->ops->destroy) {
148:     (*sp->ops->destroy)(sp);
149:     sp->ops->destroy = NULL;
150:   }
151:   (*r)(sp);
152:   PetscObjectChangeTypeName((PetscObject) sp, name);
153:   return(0);
154: }

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

159:   Not Collective

161:   Input Parameter:
162: . sp  - The PetscDualSpace

164:   Output Parameter:
165: . name - The PetscDualSpace type name

167:   Level: intermediate

169: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
170: @*/
171: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
172: {

178:   if (!PetscDualSpaceRegisterAllCalled) {
179:     PetscDualSpaceRegisterAll();
180:   }
181:   *name = ((PetscObject) sp)->type_name;
182:   return(0);
183: }

185: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
186: {
187:   PetscViewerFormat format;
188:   PetscInt          pdim, f;
189:   PetscErrorCode    ierr;

192:   PetscDualSpaceGetDimension(sp, &pdim);
193:   PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);
194:   PetscViewerASCIIPushTab(v);
195:   PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
196:   if (sp->ops->view) {(*sp->ops->view)(sp, v);}
197:   PetscViewerGetFormat(v, &format);
198:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
199:     PetscViewerASCIIPushTab(v);
200:     for (f = 0; f < pdim; ++f) {
201:       PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
202:       PetscViewerASCIIPushTab(v);
203:       PetscQuadratureView(sp->functional[f], v);
204:       PetscViewerASCIIPopTab(v);
205:     }
206:     PetscViewerASCIIPopTab(v);
207:   }
208:   PetscViewerASCIIPopTab(v);
209:   return(0);
210: }

212: /*@
213:   PetscDualSpaceView - Views a PetscDualSpace

215:   Collective on sp

217:   Input Parameter:
218: + sp - the PetscDualSpace object to view
219: - v  - the viewer

221:   Level: developer

223: .seealso PetscDualSpaceDestroy()
224: @*/
225: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
226: {
227:   PetscBool      iascii;

233:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
234:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
235:   if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
236:   return(0);
237: }

239: /*@
240:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

242:   Collective on sp

244:   Input Parameter:
245: . sp - the PetscDualSpace object to set options for

247:   Options Database:
248: . -petscspace_degree the approximation order of the space

250:   Level: developer

252: .seealso PetscDualSpaceView()
253: @*/
254: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
255: {
256:   PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
257:   PetscInt                    refDim  = 0;
258:   PetscBool                   flg;
259:   const char                 *defaultType;
260:   char                        name[256];
261:   PetscErrorCode              ierr;

265:   if (!((PetscObject) sp)->type_name) {
266:     defaultType = PETSCDUALSPACELAGRANGE;
267:   } else {
268:     defaultType = ((PetscObject) sp)->type_name;
269:   }
270:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

272:   PetscObjectOptionsBegin((PetscObject) sp);
273:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
274:   if (flg) {
275:     PetscDualSpaceSetType(sp, name);
276:   } else if (!((PetscObject) sp)->type_name) {
277:     PetscDualSpaceSetType(sp, defaultType);
278:   }
279:   PetscOptionsInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
280:   PetscOptionsInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL);
281:   if (sp->ops->setfromoptions) {
282:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
283:   }
284:   PetscOptionsInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL);
285:   PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
286:   if (flg) {
287:     DM K;

289:     if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
290:     PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
291:     PetscDualSpaceSetDM(sp, K);
292:     DMDestroy(&K);
293:   }

295:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
296:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
297:   PetscOptionsEnd();
298:   sp->setfromoptionscalled = PETSC_TRUE;
299:   return(0);
300: }

302: /*@
303:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

305:   Collective on sp

307:   Input Parameter:
308: . sp - the PetscDualSpace object to setup

310:   Level: developer

312: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
313: @*/
314: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
315: {

320:   if (sp->setupcalled) return(0);
321:   sp->setupcalled = PETSC_TRUE;
322:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
323:   if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
324:   return(0);
325: }

327: /*@
328:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

330:   Collective on sp

332:   Input Parameter:
333: . sp - the PetscDualSpace object to destroy

335:   Level: developer

337: .seealso PetscDualSpaceView()
338: @*/
339: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
340: {
341:   PetscInt       dim, f;

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

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

351:   PetscDualSpaceGetDimension(*sp, &dim);
352:   for (f = 0; f < dim; ++f) {
353:     PetscQuadratureDestroy(&(*sp)->functional[f]);
354:   }
355:   PetscFree((*sp)->functional);
356:   PetscQuadratureDestroy(&(*sp)->allPoints);
357:   DMDestroy(&(*sp)->dm);

359:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
360:   PetscHeaderDestroy(sp);
361:   return(0);
362: }

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

367:   Collective

369:   Input Parameter:
370: . comm - The communicator for the PetscDualSpace object

372:   Output Parameter:
373: . sp - The PetscDualSpace object

375:   Level: beginner

377: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
378: @*/
379: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
380: {
381:   PetscDualSpace s;

386:   PetscCitationsRegister(FECitation,&FEcite);
387:   *sp  = NULL;
388:   PetscFEInitializePackage();

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

392:   s->order = 0;
393:   s->Nc    = 1;
394:   s->k     = 0;
395:   s->setupcalled = PETSC_FALSE;

397:   *sp = s;
398:   return(0);
399: }

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

404:   Collective on sp

406:   Input Parameter:
407: . sp - The original PetscDualSpace

409:   Output Parameter:
410: . spNew - The duplicate PetscDualSpace

412:   Level: beginner

414: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
415: @*/
416: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
417: {

423:   (*sp->ops->duplicate)(sp, spNew);
424:   return(0);
425: }

427: /*@
428:   PetscDualSpaceGetDM - Get the DM representing the reference cell

430:   Not collective

432:   Input Parameter:
433: . sp - The PetscDualSpace

435:   Output Parameter:
436: . dm - The reference cell

438:   Level: intermediate

440: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
441: @*/
442: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
443: {
447:   *dm = sp->dm;
448:   return(0);
449: }

451: /*@
452:   PetscDualSpaceSetDM - Get the DM representing the reference cell

454:   Not collective

456:   Input Parameters:
457: + sp - The PetscDualSpace
458: - dm - The reference cell

460:   Level: intermediate

462: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
463: @*/
464: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
465: {

471:   DMDestroy(&sp->dm);
472:   PetscObjectReference((PetscObject) dm);
473:   sp->dm = dm;
474:   return(0);
475: }

477: /*@
478:   PetscDualSpaceGetOrder - Get the order of the dual space

480:   Not collective

482:   Input Parameter:
483: . sp - The PetscDualSpace

485:   Output Parameter:
486: . order - The order

488:   Level: intermediate

490: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
491: @*/
492: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
493: {
497:   *order = sp->order;
498:   return(0);
499: }

501: /*@
502:   PetscDualSpaceSetOrder - Set the order of the dual space

504:   Not collective

506:   Input Parameters:
507: + sp - The PetscDualSpace
508: - order - The order

510:   Level: intermediate

512: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
513: @*/
514: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
515: {
518:   sp->order = order;
519:   return(0);
520: }

522: /*@
523:   PetscDualSpaceGetNumComponents - Return the number of components for this space

525:   Input Parameter:
526: . sp - The PetscDualSpace

528:   Output Parameter:
529: . Nc - The number of components

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

533:   Level: intermediate

535: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
536: @*/
537: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
538: {
542:   *Nc = sp->Nc;
543:   return(0);
544: }

546: /*@
547:   PetscDualSpaceSetNumComponents - Set the number of components for this space

549:   Input Parameters:
550: + sp - The PetscDualSpace
551: - order - The number of components

553:   Level: intermediate

555: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
556: @*/
557: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
558: {
561:   sp->Nc = Nc;
562:   return(0);
563: }

565: /*@
566:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

568:   Not collective

570:   Input Parameters:
571: + sp - The PetscDualSpace
572: - i  - The basis number

574:   Output Parameter:
575: . functional - The basis functional

577:   Level: intermediate

579: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
580: @*/
581: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
582: {
583:   PetscInt       dim;

589:   PetscDualSpaceGetDimension(sp, &dim);
590:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
591:   *functional = sp->functional[i];
592:   return(0);
593: }

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

598:   Not collective

600:   Input Parameter:
601: . sp - The PetscDualSpace

603:   Output Parameter:
604: . dim - The dimension

606:   Level: intermediate

608: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
609: @*/
610: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
611: {

617:   *dim = 0;
618:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
619:   return(0);
620: }

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

625:   Not collective

627:   Input Parameter:
628: . sp - The PetscDualSpace

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

633:   Level: intermediate

635: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
636: @*/
637: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
638: {

644:   (*sp->ops->getnumdof)(sp, numDof);
645:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
646:   return(0);
647: }

649: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
650: {
651:   DM             dm;
652:   PetscInt       pStart, pEnd, depth, h, offset;
653:   const PetscInt *numDof;

657:   PetscDualSpaceGetDM(sp,&dm);
658:   DMPlexGetChart(dm,&pStart,&pEnd);
659:   PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
660:   PetscSectionSetChart(*section,pStart,pEnd);
661:   DMPlexGetDepth(dm,&depth);
662:   PetscDualSpaceGetNumDof(sp,&numDof);
663:   for (h = 0; h <= depth; h++) {
664:     PetscInt hStart, hEnd, p, dof;

666:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
667:     dof = numDof[depth - h];
668:     for (p = hStart; p < hEnd; p++) {
669:       PetscSectionSetDof(*section,p,dof);
670:     }
671:   }
672:   PetscSectionSetUp(*section);
673:   for (h = 0, offset = 0; h <= depth; h++) {
674:     PetscInt hStart, hEnd, p, dof;

676:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
677:     dof = numDof[depth - h];
678:     for (p = hStart; p < hEnd; p++) {
679:       PetscSectionGetDof(*section,p,&dof);
680:       PetscSectionSetOffset(*section,p,offset);
681:       offset += dof;
682:     }
683:   }
684:   return(0);
685: }

687: /*@
688:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

690:   Collective on sp

692:   Input Parameters:
693: + sp      - The PetscDualSpace
694: . dim     - The spatial dimension
695: - simplex - Flag for simplex, otherwise use a tensor-product cell

697:   Output Parameter:
698: . refdm - The reference cell

700:   Level: advanced

702: .seealso: PetscDualSpaceCreate(), DMPLEX
703: @*/
704: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
705: {

709:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
710:   return(0);
711: }

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

716:   Input Parameters:
717: + sp      - The PetscDualSpace object
718: . f       - The basis functional index
719: . time    - The time
720: . 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)
721: . numComp - The number of components for the function
722: . func    - The input function
723: - ctx     - A context for the function

725:   Output Parameter:
726: . value   - numComp output values

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

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

733:   Level: developer

735: .seealso: PetscDualSpaceCreate()
736: @*/
737: 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)
738: {

745:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
746:   return(0);
747: }

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

752:   Input Parameters:
753: + sp        - The PetscDualSpace object
754: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

756:   Output Parameter:
757: . spValue   - The values of all dual space functionals

759:   Level: developer

761: .seealso: PetscDualSpaceCreate()
762: @*/
763: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
764: {

769:   (*sp->ops->applyall)(sp, pointEval, spValue);
770:   return(0);
771: }

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

776:   Input Parameters:
777: + sp    - The PetscDualSpace object
778: . f     - The basis functional index
779: . time  - The time
780: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
781: . Nc    - The number of components for the function
782: . func  - The input function
783: - ctx   - A context for the function

785:   Output Parameter:
786: . value   - The output value

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

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

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

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

797: where both n and f have Nc components.

799:   Level: developer

801: .seealso: PetscDualSpaceCreate()
802: @*/
803: 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)
804: {
805:   DM               dm;
806:   PetscQuadrature  n;
807:   const PetscReal *points, *weights;
808:   PetscReal        x[3];
809:   PetscScalar     *val;
810:   PetscInt         dim, dE, qNc, c, Nq, q;
811:   PetscBool        isAffine;
812:   PetscErrorCode   ierr;

817:   PetscDualSpaceGetDM(sp, &dm);
818:   PetscDualSpaceGetFunctional(sp, f, &n);
819:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
820:   if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
821:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
822:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
823:   *value = 0.0;
824:   isAffine = cgeom->isAffine;
825:   dE = cgeom->dimEmbed;
826:   for (q = 0; q < Nq; ++q) {
827:     if (isAffine) {
828:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
829:       (*func)(dE, time, x, Nc, val, ctx);
830:     } else {
831:       (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
832:     }
833:     for (c = 0; c < Nc; ++c) {
834:       *value += val[c]*weights[q*Nc+c];
835:     }
836:   }
837:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
838:   return(0);
839: }

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

844:   Input Parameters:
845: + sp        - The PetscDualSpace object
846: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

848:   Output Parameter:
849: . spValue   - The values of all dual space functionals

851:   Level: developer

853: .seealso: PetscDualSpaceCreate()
854: @*/
855: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
856: {
857:   PetscQuadrature  n;
858:   const PetscReal *points, *weights;
859:   PetscInt         qNc, c, Nq, q, f, spdim, Nc;
860:   PetscInt         offset;
861:   PetscErrorCode   ierr;

867:   PetscDualSpaceGetDimension(sp, &spdim);
868:   PetscDualSpaceGetNumComponents(sp, &Nc);
869:   for (f = 0, offset = 0; f < spdim; f++) {
870:     PetscDualSpaceGetFunctional(sp, f, &n);
871:     PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
872:     if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
873:     spValue[f] = 0.0;
874:     for (q = 0; q < Nq; ++q) {
875:       for (c = 0; c < Nc; ++c) {
876:         spValue[f] += pointEval[offset++]*weights[q*Nc+c];
877:       }
878:     }
879:   }
880:   return(0);
881: }

883: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
884: {

890:   if (!sp->allPoints && sp->ops->createallpoints) {
891:     (*sp->ops->createallpoints)(sp,&sp->allPoints);
892:   }
893:   *allPoints = sp->allPoints;
894:   return(0);
895: }

897: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
898: {
899:   PetscInt        spdim;
900:   PetscInt        numPoints, offset;
901:   PetscReal       *points;
902:   PetscInt        f, dim;
903:   PetscQuadrature q;
904:   PetscErrorCode  ierr;

907:   PetscDualSpaceGetDimension(sp,&spdim);
908:   if (!spdim) {
909:     PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
910:     PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
911:   }
912:   PetscDualSpaceGetFunctional(sp,0,&q);
913:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
914:   for (f = 1; f < spdim; f++) {
915:     PetscInt Np;

917:     PetscDualSpaceGetFunctional(sp,f,&q);
918:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
919:     numPoints += Np;
920:   }
921:   PetscMalloc1(dim*numPoints,&points);
922:   for (f = 0, offset = 0; f < spdim; f++) {
923:     const PetscReal *p;
924:     PetscInt        Np, i;

926:     PetscDualSpaceGetFunctional(sp,f,&q);
927:     PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
928:     for (i = 0; i < Np * dim; i++) {
929:       points[offset + i] = p[i];
930:     }
931:     offset += Np * dim;
932:   }
933:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
934:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
935:   return(0);
936: }

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

941:   Input Parameters:
942: + sp    - The PetscDualSpace object
943: . f     - The basis functional index
944: . time  - The time
945: . cgeom - A context with geometric information for this cell, we currently just use the centroid
946: . Nc    - The number of components for the function
947: . func  - The input function
948: - ctx   - A context for the function

950:   Output Parameter:
951: . value - The output value (scalar)

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

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

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

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

962: where both n and f have Nc components.

964:   Level: developer

966: .seealso: PetscDualSpaceCreate()
967: @*/
968: 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)
969: {
970:   DM               dm;
971:   PetscQuadrature  n;
972:   const PetscReal *points, *weights;
973:   PetscScalar     *val;
974:   PetscInt         dimEmbed, qNc, c, Nq, q;
975:   PetscErrorCode   ierr;

980:   PetscDualSpaceGetDM(sp, &dm);
981:   DMGetCoordinateDim(dm, &dimEmbed);
982:   PetscDualSpaceGetFunctional(sp, f, &n);
983:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
984:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
985:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
986:   *value = 0.;
987:   for (q = 0; q < Nq; ++q) {
988:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
989:     for (c = 0; c < Nc; ++c) {
990:       *value += val[c]*weights[q*Nc+c];
991:     }
992:   }
993:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
994:   return(0);
995: }

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

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

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

1007:   Not collective

1009:   Input Parameters:
1010: + sp - the PetscDualSpace object
1011: - height - the height of the mesh point for which the subspace is desired

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

1017:   Level: advanced

1019: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1020: @*/
1021: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1022: {

1028:   *subsp = NULL;
1029:   if (sp->ops->getheightsubspace) {(*sp->ops->getheightsubspace)(sp, height, subsp);}
1030:   return(0);
1031: }

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

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

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

1042:   Not collective

1044:   Input Parameters:
1045: + sp - the PetscDualSpace object
1046: - point - the point (in the dual space's DM) for which the subspace is desired

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

1052:   Level: advanced

1054: .seealso: PetscDualSpace
1055: @*/
1056: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1057: {

1063:   *bdsp = NULL;
1064:   if (sp->ops->getpointsubspace) {
1065:     (*sp->ops->getpointsubspace)(sp,point,bdsp);
1066:   } else if (sp->ops->getheightsubspace) {
1067:     DM       dm;
1068:     DMLabel  label;
1069:     PetscInt dim, depth, height;

1071:     PetscDualSpaceGetDM(sp,&dm);
1072:     DMPlexGetDepth(dm,&dim);
1073:     DMPlexGetDepthLabel(dm,&label);
1074:     DMLabelGetValue(label,point,&depth);
1075:     height = dim - depth;
1076:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1077:   }
1078:   return(0);
1079: }

1081: /*@C
1082:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1084:   Not collective

1086:   Input Parameter:
1087: . sp - the PetscDualSpace object

1089:   Output Parameters:
1090: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
1091: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation

1093:   Note: The permutation and flip arrays are organized in the following way
1094: $ perms[p][ornt][dof # on point] = new local dof #
1095: $ flips[p][ornt][dof # on point] = reversal or not

1097:   Level: developer

1099: .seealso: PetscDualSpaceSetSymmetries()
1100: @*/
1101: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1102: {

1109:   if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1110:   return(0);
1111: }

1113: /*@
1114:   PetscDualSpaceGetDeRahm - Get the k-simplex associated with the functionals in this dual space

1116:   Input Parameter:
1117: . dsp - The PetscDualSpace

1119:   Output Parameter:
1120: . k   - The simplex dimension

1122:   Level: advanced

1124:   Note: Currently supported values are
1125: $ 0: These are H_1 methods that only transform coordinates
1126: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1127: $ 2: These are the same as 1
1128: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)

1130: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1131: @*/
1132: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1133: {
1137:   *k = dsp->k;
1138:   return(0);
1139: }

1141: /*@C
1142:   PetscDualSpaceTransform - Transform the function values

1144:   Input Parameters:
1145: + dsp       - The PetscDualSpace
1146: . trans     - The type of transform
1147: . isInverse - Flag to invert the transform
1148: . fegeom    - The cell geometry
1149: . Nv        - The number of function samples
1150: . Nc        - The number of function components
1151: - vals      - The function values

1153:   Output Parameter:
1154: . vals      - The transformed function values

1156:   Level: developer

1158: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1159: @*/
1160: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1161: {
1162:   PetscInt dim, v, c;

1168:   dim = dsp->dm->dim;
1169:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1170:   if (Nc == 1 || Nc != dim) return(0);
1171:   switch (trans) {
1172:     case IDENTITY_TRANSFORM: break;
1173:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1174:     if (isInverse) {
1175:       for (v = 0; v < Nv; ++v) {
1176:         switch (dim)
1177:         {
1178:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1179:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1180:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1181:         }
1182:       }
1183:     } else {
1184:       for (v = 0; v < Nv; ++v) {
1185:         switch (dim)
1186:         {
1187:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1188:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1189:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1190:         }
1191:       }
1192:     }
1193:     break;
1194:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1195:     if (isInverse) {
1196:       for (v = 0; v < Nv; ++v) {
1197:         switch (dim)
1198:         {
1199:           case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1200:           case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1201:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1202:         }
1203:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] *= fegeom->detJ[0];
1204:       }
1205:     } else {
1206:       for (v = 0; v < Nv; ++v) {
1207:         switch (dim)
1208:         {
1209:           case 2: DMPlex_Mult2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1210:           case 3: DMPlex_Mult3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1211:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1212:         }
1213:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] /= fegeom->detJ[0];
1214:       }
1215:     }
1216:     break;
1217:   }
1218:   return(0);
1219: }
1220: /*@C
1221:   PetscDualSpaceTransformGradient - Transform the function gradient values

1223:   Input Parameters:
1224: + dsp       - The PetscDualSpace
1225: . trans     - The type of transform
1226: . isInverse - Flag to invert the transform
1227: . fegeom    - The cell geometry
1228: . Nv        - The number of function gradient samples
1229: . Nc        - The number of function components
1230: - vals      - The function gradient values

1232:   Output Parameter:
1233: . vals      - The transformed function values

1235:   Level: developer

1237: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1238: @*/
1239: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1240: {
1241:   PetscInt dim, v, c, d;

1247:   dim = dsp->dm->dim;
1248:   /* Transform gradient */
1249:   for (v = 0; v < Nv; ++v) {
1250:     for (c = 0; c < Nc; ++c) {
1251:       switch (dim)
1252:       {
1253:         case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];
1254:         case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1255:         case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1256:         default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1257:       }
1258:     }
1259:   }
1260:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1261:   if (Nc == 1 || Nc != dim) return(0);
1262:   switch (trans) {
1263:     case IDENTITY_TRANSFORM: break;
1264:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1265:     if (isInverse) {
1266:       for (v = 0; v < Nv; ++v) {
1267:         for (d = 0; d < dim; ++d) {
1268:           switch (dim)
1269:           {
1270:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1271:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1272:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1273:           }
1274:         }
1275:       }
1276:     } else {
1277:       for (v = 0; v < Nv; ++v) {
1278:         for (d = 0; d < dim; ++d) {
1279:           switch (dim)
1280:           {
1281:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1282:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1283:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1284:           }
1285:         }
1286:       }
1287:     }
1288:     break;
1289:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1290:     if (isInverse) {
1291:       for (v = 0; v < Nv; ++v) {
1292:         for (d = 0; d < dim; ++d) {
1293:           switch (dim)
1294:           {
1295:             case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1296:             case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1297:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1298:           }
1299:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1300:         }
1301:       }
1302:     } else {
1303:       for (v = 0; v < Nv; ++v) {
1304:         for (d = 0; d < dim; ++d) {
1305:           switch (dim)
1306:           {
1307:             case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1308:             case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1309:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1310:           }
1311:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1312:         }
1313:       }
1314:     }
1315:     break;
1316:   }
1317:   return(0);
1318: }

1320: /*@C
1321:   PetscDualSpacePullback - Transform the given functional so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.

1323:   Input Parameters:
1324: + dsp        - The PetscDualSpace
1325: . fegeom     - The geometry for this cell
1326: . Nq         - The number of function samples
1327: . Nc         - The number of function components
1328: - pointEval  - The function values

1330:   Output Parameter:
1331: . pointEval  - The transformed function values

1333:   Level: advanced

1335:   Note: Functions transform in a complementary way (pushforward) to functionals, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.

1337: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1338: @*/
1339: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1340: {
1341:   PetscDualSpaceTransformType trans;
1342:   PetscErrorCode              ierr;

1348:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1349:      This determines their transformation properties. */
1350:   switch (dsp->k)
1351:   {
1352:     case 0: /* H^1 point evaluations */
1353:     trans = IDENTITY_TRANSFORM;break;
1354:     case 1: /* Hcurl preserves tangential edge traces  */
1355:     case 2:
1356:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1357:     case 3: /* Hdiv preserve normal traces */
1358:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1359:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1360:   }
1361:   PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
1362:   return(0);
1363: }

1365: /*@C
1366:   PetscDualSpacePushforward - Transform the given function so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.

1368:   Input Parameters:
1369: + dsp        - The PetscDualSpace
1370: . fegeom     - The geometry for this cell
1371: . Nq         - The number of function samples
1372: . Nc         - The number of function components
1373: - pointEval  - The function values

1375:   Output Parameter:
1376: . pointEval  - The transformed function values

1378:   Level: advanced

1380:   Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.

1382: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1383: @*/
1384: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1385: {
1386:   PetscDualSpaceTransformType trans;
1387:   PetscErrorCode              ierr;

1393:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1394:      This determines their transformation properties. */
1395:   switch (dsp->k)
1396:   {
1397:     case 0: /* H^1 point evaluations */
1398:     trans = IDENTITY_TRANSFORM;break;
1399:     case 1: /* Hcurl preserves tangential edge traces  */
1400:     case 2:
1401:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1402:     case 3: /* Hdiv preserve normal traces */
1403:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1404:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1405:   }
1406:   PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1407:   return(0);
1408: }

1410: /*@C
1411:   PetscDualSpacePushforwardGradient - Transform the given function gradient so that it operates on real space, rather than the reference element. Operationally, this means that we map the function evaluations depending on continuity requirements of our finite element method.

1413:   Input Parameters:
1414: + dsp        - The PetscDualSpace
1415: . fegeom     - The geometry for this cell
1416: . Nq         - The number of function gradient samples
1417: . Nc         - The number of function components
1418: - pointEval  - The function gradient values

1420:   Output Parameter:
1421: . pointEval  - The transformed function gradient values

1423:   Level: advanced

1425:   Note: Functionals transform in a complementary way (pullback) to functions, so that the scalar product is invariant. The type of transform is dependent on the associated k-simplex from the DeRahm complex.

1427: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1428: @*/PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1429: {
1430:   PetscDualSpaceTransformType trans;
1431:   PetscErrorCode              ierr;

1437:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1438:      This determines their transformation properties. */
1439:   switch (dsp->k)
1440:   {
1441:     case 0: /* H^1 point evaluations */
1442:     trans = IDENTITY_TRANSFORM;break;
1443:     case 1: /* Hcurl preserves tangential edge traces  */
1444:     case 2:
1445:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1446:     case 3: /* Hdiv preserve normal traces */
1447:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1448:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1449:   }
1450:   PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1451:   return(0);
1452: }