Actual source code: dualspace.c

petsc-master 2019-12-13
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: /*@C
213:    PetscDualSpaceViewFromOptions - View from Options

215:    Collective on PetscDualSpace

217:    Input Parameters:
218: +  A - the PetscDualSpace object
219: .  obj - Optional object, proivides prefix
220: -  name - command line option

222:    Level: intermediate
223: .seealso:  PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate()
224: @*/
225: PetscErrorCode  PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[])
226: {

231:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
232:   return(0);
233: }

235: /*@
236:   PetscDualSpaceView - Views a PetscDualSpace

238:   Collective on sp

240:   Input Parameter:
241: + sp - the PetscDualSpace object to view
242: - v  - the viewer

244:   Level: beginner

246: .seealso PetscDualSpaceDestroy(), PetscDualSpace
247: @*/
248: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
249: {
250:   PetscBool      iascii;

256:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
257:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
258:   if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
259:   return(0);
260: }

262: /*@
263:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

265:   Collective on sp

267:   Input Parameter:
268: . sp - the PetscDualSpace object to set options for

270:   Options Database:
271: . -petscspace_degree the approximation order of the space

273:   Level: intermediate

275: .seealso PetscDualSpaceView(), PetscDualSpace, PetscObjectSetFromOptions()
276: @*/
277: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
278: {
279:   PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
280:   PetscInt                    refDim  = 0;
281:   PetscBool                   flg;
282:   const char                 *defaultType;
283:   char                        name[256];
284:   PetscErrorCode              ierr;

288:   if (!((PetscObject) sp)->type_name) {
289:     defaultType = PETSCDUALSPACELAGRANGE;
290:   } else {
291:     defaultType = ((PetscObject) sp)->type_name;
292:   }
293:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

295:   PetscObjectOptionsBegin((PetscObject) sp);
296:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
297:   if (flg) {
298:     PetscDualSpaceSetType(sp, name);
299:   } else if (!((PetscObject) sp)->type_name) {
300:     PetscDualSpaceSetType(sp, defaultType);
301:   }
302:   PetscOptionsBoundedInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
303:   PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
304:   if (sp->ops->setfromoptions) {
305:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
306:   }
307:   PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
308:   PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
309:   if (flg) {
310:     DM K;

312:     if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
313:     PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
314:     PetscDualSpaceSetDM(sp, K);
315:     DMDestroy(&K);
316:   }

318:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
319:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
320:   PetscOptionsEnd();
321:   sp->setfromoptionscalled = PETSC_TRUE;
322:   return(0);
323: }

325: /*@
326:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

328:   Collective on sp

330:   Input Parameter:
331: . sp - the PetscDualSpace object to setup

333:   Level: intermediate

335: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace
336: @*/
337: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
338: {

343:   if (sp->setupcalled) return(0);
344:   sp->setupcalled = PETSC_TRUE;
345:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
346:   if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
347:   return(0);
348: }

350: /*@
351:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

353:   Collective on sp

355:   Input Parameter:
356: . sp - the PetscDualSpace object to destroy

358:   Level: beginner

360: .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate()
361: @*/
362: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
363: {
364:   PetscInt       dim, f;

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

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

374:   PetscDualSpaceGetDimension(*sp, &dim);
375:   for (f = 0; f < dim; ++f) {
376:     PetscQuadratureDestroy(&(*sp)->functional[f]);
377:   }
378:   PetscFree((*sp)->functional);
379:   PetscQuadratureDestroy(&(*sp)->allPoints);
380:   DMDestroy(&(*sp)->dm);

382:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
383:   PetscHeaderDestroy(sp);
384:   return(0);
385: }

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

390:   Collective

392:   Input Parameter:
393: . comm - The communicator for the PetscDualSpace object

395:   Output Parameter:
396: . sp - The PetscDualSpace object

398:   Level: beginner

400: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
401: @*/
402: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
403: {
404:   PetscDualSpace s;

409:   PetscCitationsRegister(FECitation,&FEcite);
410:   *sp  = NULL;
411:   PetscFEInitializePackage();

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

415:   s->order = 0;
416:   s->Nc    = 1;
417:   s->k     = 0;
418:   s->setupcalled = PETSC_FALSE;

420:   *sp = s;
421:   return(0);
422: }

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

427:   Collective on sp

429:   Input Parameter:
430: . sp - The original PetscDualSpace

432:   Output Parameter:
433: . spNew - The duplicate PetscDualSpace

435:   Level: beginner

437: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
438: @*/
439: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
440: {

446:   (*sp->ops->duplicate)(sp, spNew);
447:   return(0);
448: }

450: /*@
451:   PetscDualSpaceGetDM - Get the DM representing the reference cell

453:   Not collective

455:   Input Parameter:
456: . sp - The PetscDualSpace

458:   Output Parameter:
459: . dm - The reference cell

461:   Level: intermediate

463: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
464: @*/
465: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
466: {
470:   *dm = sp->dm;
471:   return(0);
472: }

474: /*@
475:   PetscDualSpaceSetDM - Get the DM representing the reference cell

477:   Not collective

479:   Input Parameters:
480: + sp - The PetscDualSpace
481: - dm - The reference cell

483:   Level: intermediate

485: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
486: @*/
487: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
488: {

494:   DMDestroy(&sp->dm);
495:   PetscObjectReference((PetscObject) dm);
496:   sp->dm = dm;
497:   return(0);
498: }

500: /*@
501:   PetscDualSpaceGetOrder - Get the order of the dual space

503:   Not collective

505:   Input Parameter:
506: . sp - The PetscDualSpace

508:   Output Parameter:
509: . order - The order

511:   Level: intermediate

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

524: /*@
525:   PetscDualSpaceSetOrder - Set the order of the dual space

527:   Not collective

529:   Input Parameters:
530: + sp - The PetscDualSpace
531: - order - The order

533:   Level: intermediate

535: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
536: @*/
537: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
538: {
541:   sp->order = order;
542:   return(0);
543: }

545: /*@
546:   PetscDualSpaceGetNumComponents - Return the number of components for this space

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

551:   Output Parameter:
552: . Nc - The number of components

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

556:   Level: intermediate

558: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
559: @*/
560: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
561: {
565:   *Nc = sp->Nc;
566:   return(0);
567: }

569: /*@
570:   PetscDualSpaceSetNumComponents - Set the number of components for this space

572:   Input Parameters:
573: + sp - The PetscDualSpace
574: - order - The number of components

576:   Level: intermediate

578: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
579: @*/
580: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
581: {
584:   sp->Nc = Nc;
585:   return(0);
586: }

588: /*@
589:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

591:   Not collective

593:   Input Parameters:
594: + sp - The PetscDualSpace
595: - i  - The basis number

597:   Output Parameter:
598: . functional - The basis functional

600:   Level: intermediate

602: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
603: @*/
604: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
605: {
606:   PetscInt       dim;

612:   PetscDualSpaceGetDimension(sp, &dim);
613:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
614:   *functional = sp->functional[i];
615:   return(0);
616: }

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

621:   Not collective

623:   Input Parameter:
624: . sp - The PetscDualSpace

626:   Output Parameter:
627: . dim - The dimension

629:   Level: intermediate

631: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
632: @*/
633: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
634: {

640:   *dim = 0;
641:   if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
642:   return(0);
643: }

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

648:   Not collective

650:   Input Parameter:
651: . sp - The PetscDualSpace

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

656:   Level: intermediate

658: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
659: @*/
660: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
661: {

667:   (*sp->ops->getnumdof)(sp, numDof);
668:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
669:   return(0);
670: }

672: /*@
673:   PetscDualSpaceCreateSection - Create a PetscSection over the reference cell with the layout from this space

675:   Collective on sp

677:   Input Parameters:
678: + sp      - The PetscDualSpace

680:   Output Parameter:
681: . section - The section

683:   Level: advanced

685: .seealso: PetscDualSpaceCreate(), DMPLEX
686: @*/
687: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
688: {
689:   DM             dm;
690:   PetscInt       pStart, pEnd, depth, h, offset;
691:   const PetscInt *numDof;

695:   PetscDualSpaceGetDM(sp,&dm);
696:   DMPlexGetChart(dm,&pStart,&pEnd);
697:   PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
698:   PetscSectionSetChart(*section,pStart,pEnd);
699:   DMPlexGetDepth(dm,&depth);
700:   PetscDualSpaceGetNumDof(sp,&numDof);
701:   for (h = 0; h <= depth; h++) {
702:     PetscInt hStart, hEnd, p, dof;

704:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
705:     dof = numDof[depth - h];
706:     for (p = hStart; p < hEnd; p++) {
707:       PetscSectionSetDof(*section,p,dof);
708:     }
709:   }
710:   PetscSectionSetUp(*section);
711:   for (h = 0, offset = 0; h <= depth; h++) {
712:     PetscInt hStart, hEnd, p, dof;

714:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
715:     dof = numDof[depth - h];
716:     for (p = hStart; p < hEnd; p++) {
717:       PetscSectionGetDof(*section,p,&dof);
718:       PetscSectionSetOffset(*section,p,offset);
719:       offset += dof;
720:     }
721:   }
722:   return(0);
723: }

725: /*@
726:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

728:   Collective on sp

730:   Input Parameters:
731: + sp      - The PetscDualSpace
732: . dim     - The spatial dimension
733: - simplex - Flag for simplex, otherwise use a tensor-product cell

735:   Output Parameter:
736: . refdm - The reference cell

738:   Level: intermediate

740: .seealso: PetscDualSpaceCreate(), DMPLEX
741: @*/
742: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
743: {

747:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
748:   return(0);
749: }

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

754:   Input Parameters:
755: + sp      - The PetscDualSpace object
756: . f       - The basis functional index
757: . time    - The time
758: . 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)
759: . numComp - The number of components for the function
760: . func    - The input function
761: - ctx     - A context for the function

763:   Output Parameter:
764: . value   - numComp output values

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

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

771:   Level: beginner

773: .seealso: PetscDualSpaceCreate()
774: @*/
775: 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)
776: {

783:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
784:   return(0);
785: }

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

790:   Input Parameters:
791: + sp        - The PetscDualSpace object
792: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

794:   Output Parameter:
795: . spValue   - The values of all dual space functionals

797:   Level: beginner

799: .seealso: PetscDualSpaceCreate()
800: @*/
801: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
802: {

807:   (*sp->ops->applyall)(sp, pointEval, spValue);
808:   return(0);
809: }

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

814:   Input Parameters:
815: + sp    - The PetscDualSpace object
816: . f     - The basis functional index
817: . time  - The time
818: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
819: . Nc    - The number of components for the function
820: . func  - The input function
821: - ctx   - A context for the function

823:   Output Parameter:
824: . value   - The output value

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

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

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

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

835: where both n and f have Nc components.

837:   Level: beginner

839: .seealso: PetscDualSpaceCreate()
840: @*/
841: 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)
842: {
843:   DM               dm;
844:   PetscQuadrature  n;
845:   const PetscReal *points, *weights;
846:   PetscReal        x[3];
847:   PetscScalar     *val;
848:   PetscInt         dim, dE, qNc, c, Nq, q;
849:   PetscBool        isAffine;
850:   PetscErrorCode   ierr;

855:   PetscDualSpaceGetDM(sp, &dm);
856:   PetscDualSpaceGetFunctional(sp, f, &n);
857:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
858:   if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
859:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
860:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
861:   *value = 0.0;
862:   isAffine = cgeom->isAffine;
863:   dE = cgeom->dimEmbed;
864:   for (q = 0; q < Nq; ++q) {
865:     if (isAffine) {
866:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
867:       (*func)(dE, time, x, Nc, val, ctx);
868:     } else {
869:       (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
870:     }
871:     for (c = 0; c < Nc; ++c) {
872:       *value += val[c]*weights[q*Nc+c];
873:     }
874:   }
875:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
876:   return(0);
877: }

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

882:   Input Parameters:
883: + sp        - The PetscDualSpace object
884: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

886:   Output Parameter:
887: . spValue   - The values of all dual space functionals

889:   Level: beginner

891: .seealso: PetscDualSpaceCreate()
892: @*/
893: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
894: {
895:   PetscQuadrature  n;
896:   const PetscReal *points, *weights;
897:   PetscInt         qNc, c, Nq, q, f, spdim, Nc;
898:   PetscInt         offset;
899:   PetscErrorCode   ierr;

905:   PetscDualSpaceGetDimension(sp, &spdim);
906:   PetscDualSpaceGetNumComponents(sp, &Nc);
907:   for (f = 0, offset = 0; f < spdim; f++) {
908:     PetscDualSpaceGetFunctional(sp, f, &n);
909:     PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
910:     if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
911:     spValue[f] = 0.0;
912:     for (q = 0; q < Nq; ++q) {
913:       for (c = 0; c < Nc; ++c) {
914:         spValue[f] += pointEval[offset++]*weights[q*Nc+c];
915:       }
916:     }
917:   }
918:   return(0);
919: }

921: /*@
922:   PetscDualSpaceGetAllPoints - Get all quadrature points from this space

924:   Input Parameter:
925: . sp - The dualspace

927:   Output Parameter:
928: . allPoints - A PetscQuadrature object containing all evaluation points

930:   Level: advanced

932: .seealso: PetscDualSpaceCreate()
933: @*/
934: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
935: {

941:   if (!sp->allPoints && sp->ops->createallpoints) {
942:     (*sp->ops->createallpoints)(sp,&sp->allPoints);
943:   }
944:   *allPoints = sp->allPoints;
945:   return(0);
946: }

948: /*@
949:   PetscDualSpaceCreateAllPointsDefault - Create all evaluation points by examining functionals

951:   Input Parameter:
952: . sp - The dualspace

954:   Output Parameter:
955: . allPoints - A PetscQuadrature object containing all evaluation points

957:   Level: advanced

959: .seealso: PetscDualSpaceCreate()
960: @*/
961: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
962: {
963:   PetscInt        spdim;
964:   PetscInt        numPoints, offset;
965:   PetscReal       *points;
966:   PetscInt        f, dim;
967:   PetscQuadrature q;
968:   PetscErrorCode  ierr;

971:   PetscDualSpaceGetDimension(sp,&spdim);
972:   if (!spdim) {
973:     PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
974:     PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
975:   }
976:   PetscDualSpaceGetFunctional(sp,0,&q);
977:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
978:   for (f = 1; f < spdim; f++) {
979:     PetscInt Np;

981:     PetscDualSpaceGetFunctional(sp,f,&q);
982:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
983:     numPoints += Np;
984:   }
985:   PetscMalloc1(dim*numPoints,&points);
986:   for (f = 0, offset = 0; f < spdim; f++) {
987:     const PetscReal *p;
988:     PetscInt        Np, i;

990:     PetscDualSpaceGetFunctional(sp,f,&q);
991:     PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
992:     for (i = 0; i < Np * dim; i++) {
993:       points[offset + i] = p[i];
994:     }
995:     offset += Np * dim;
996:   }
997:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
998:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
999:   return(0);
1000: }

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

1005:   Input Parameters:
1006: + sp    - The PetscDualSpace object
1007: . f     - The basis functional index
1008: . time  - The time
1009: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1010: . Nc    - The number of components for the function
1011: . func  - The input function
1012: - ctx   - A context for the function

1014:   Output Parameter:
1015: . value - The output value (scalar)

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

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

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

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

1026: where both n and f have Nc components.

1028:   Level: beginner

1030: .seealso: PetscDualSpaceCreate()
1031: @*/
1032: 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)
1033: {
1034:   DM               dm;
1035:   PetscQuadrature  n;
1036:   const PetscReal *points, *weights;
1037:   PetscScalar     *val;
1038:   PetscInt         dimEmbed, qNc, c, Nq, q;
1039:   PetscErrorCode   ierr;

1044:   PetscDualSpaceGetDM(sp, &dm);
1045:   DMGetCoordinateDim(dm, &dimEmbed);
1046:   PetscDualSpaceGetFunctional(sp, f, &n);
1047:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1048:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1049:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1050:   *value = 0.;
1051:   for (q = 0; q < Nq; ++q) {
1052:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1053:     for (c = 0; c < Nc; ++c) {
1054:       *value += val[c]*weights[q*Nc+c];
1055:     }
1056:   }
1057:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1058:   return(0);
1059: }

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

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

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

1071:   Not collective

1073:   Input Parameters:
1074: + sp - the PetscDualSpace object
1075: - height - the height of the mesh point for which the subspace is desired

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

1081:   Level: advanced

1083: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1084: @*/
1085: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1086: {

1092:   *subsp = NULL;
1093:   if (sp->ops->getheightsubspace) {(*sp->ops->getheightsubspace)(sp, height, subsp);}
1094:   return(0);
1095: }

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

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

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

1106:   Not collective

1108:   Input Parameters:
1109: + sp - the PetscDualSpace object
1110: - point - the point (in the dual space's DM) for which the subspace is desired

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

1116:   Level: advanced

1118: .seealso: PetscDualSpace
1119: @*/
1120: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1121: {

1127:   *bdsp = NULL;
1128:   if (sp->ops->getpointsubspace) {
1129:     (*sp->ops->getpointsubspace)(sp,point,bdsp);
1130:   } else if (sp->ops->getheightsubspace) {
1131:     DM       dm;
1132:     DMLabel  label;
1133:     PetscInt dim, depth, height;

1135:     PetscDualSpaceGetDM(sp,&dm);
1136:     DMPlexGetDepth(dm,&dim);
1137:     DMPlexGetDepthLabel(dm,&label);
1138:     DMLabelGetValue(label,point,&depth);
1139:     height = dim - depth;
1140:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1141:   }
1142:   return(0);
1143: }

1145: /*@C
1146:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1148:   Not collective

1150:   Input Parameter:
1151: . sp - the PetscDualSpace object

1153:   Output Parameters:
1154: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
1155: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation

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

1161:   Level: developer

1163: .seealso: PetscDualSpaceSetSymmetries()
1164: @*/
1165: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1166: {

1173:   if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1174:   return(0);
1175: }

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

1180:   Input Parameter:
1181: . dsp - The PetscDualSpace

1183:   Output Parameter:
1184: . k   - The simplex dimension

1186:   Level: developer

1188:   Note: Currently supported values are
1189: $ 0: These are H_1 methods that only transform coordinates
1190: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1191: $ 2: These are the same as 1
1192: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)

1194: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1195: @*/
1196: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1197: {
1201:   *k = dsp->k;
1202:   return(0);
1203: }

1205: /*@C
1206:   PetscDualSpaceTransform - Transform the function values

1208:   Input Parameters:
1209: + dsp       - The PetscDualSpace
1210: . trans     - The type of transform
1211: . isInverse - Flag to invert the transform
1212: . fegeom    - The cell geometry
1213: . Nv        - The number of function samples
1214: . Nc        - The number of function components
1215: - vals      - The function values

1217:   Output Parameter:
1218: . vals      - The transformed function values

1220:   Level: intermediate

1222: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1223: @*/
1224: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1225: {
1226:   PetscInt dim, v, c;

1232:   dim = dsp->dm->dim;
1233:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1234:   if (Nc == 1 || Nc != dim) return(0);
1235:   switch (trans) {
1236:     case IDENTITY_TRANSFORM: break;
1237:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1238:     if (isInverse) {
1239:       for (v = 0; v < Nv; ++v) {
1240:         switch (dim)
1241:         {
1242:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1243:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1244:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1245:         }
1246:       }
1247:     } else {
1248:       for (v = 0; v < Nv; ++v) {
1249:         switch (dim)
1250:         {
1251:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1252:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1253:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1254:         }
1255:       }
1256:     }
1257:     break;
1258:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1259:     if (isInverse) {
1260:       for (v = 0; v < Nv; ++v) {
1261:         switch (dim)
1262:         {
1263:           case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1264:           case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1265:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1266:         }
1267:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] *= fegeom->detJ[0];
1268:       }
1269:     } else {
1270:       for (v = 0; v < Nv; ++v) {
1271:         switch (dim)
1272:         {
1273:           case 2: DMPlex_Mult2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1274:           case 3: DMPlex_Mult3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1275:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1276:         }
1277:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] /= fegeom->detJ[0];
1278:       }
1279:     }
1280:     break;
1281:   }
1282:   return(0);
1283: }
1284: /*@C
1285:   PetscDualSpaceTransformGradient - Transform the function gradient values

1287:   Input Parameters:
1288: + dsp       - The PetscDualSpace
1289: . trans     - The type of transform
1290: . isInverse - Flag to invert the transform
1291: . fegeom    - The cell geometry
1292: . Nv        - The number of function gradient samples
1293: . Nc        - The number of function components
1294: - vals      - The function gradient values

1296:   Output Parameter:
1297: . vals      - The transformed function values

1299:   Level: intermediate

1301: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1302: @*/
1303: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1304: {
1305:   PetscInt dim, v, c, d;

1311:   dim = dsp->dm->dim;
1312:   /* Transform gradient */
1313:   for (v = 0; v < Nv; ++v) {
1314:     for (c = 0; c < Nc; ++c) {
1315:       switch (dim)
1316:       {
1317:         case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];
1318:         case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1319:         case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1320:         default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1321:       }
1322:     }
1323:   }
1324:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1325:   if (Nc == 1 || Nc != dim) return(0);
1326:   switch (trans) {
1327:     case IDENTITY_TRANSFORM: break;
1328:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1329:     if (isInverse) {
1330:       for (v = 0; v < Nv; ++v) {
1331:         for (d = 0; d < dim; ++d) {
1332:           switch (dim)
1333:           {
1334:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1335:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1336:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1337:           }
1338:         }
1339:       }
1340:     } else {
1341:       for (v = 0; v < Nv; ++v) {
1342:         for (d = 0; d < dim; ++d) {
1343:           switch (dim)
1344:           {
1345:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1346:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1347:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1348:           }
1349:         }
1350:       }
1351:     }
1352:     break;
1353:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1354:     if (isInverse) {
1355:       for (v = 0; v < Nv; ++v) {
1356:         for (d = 0; d < dim; ++d) {
1357:           switch (dim)
1358:           {
1359:             case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1360:             case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1361:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1362:           }
1363:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1364:         }
1365:       }
1366:     } else {
1367:       for (v = 0; v < Nv; ++v) {
1368:         for (d = 0; d < dim; ++d) {
1369:           switch (dim)
1370:           {
1371:             case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1372:             case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1373:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1374:           }
1375:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1376:         }
1377:       }
1378:     }
1379:     break;
1380:   }
1381:   return(0);
1382: }

1384: /*@C
1385:   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.

1387:   Input Parameters:
1388: + dsp        - The PetscDualSpace
1389: . fegeom     - The geometry for this cell
1390: . Nq         - The number of function samples
1391: . Nc         - The number of function components
1392: - pointEval  - The function values

1394:   Output Parameter:
1395: . pointEval  - The transformed function values

1397:   Level: advanced

1399:   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.

1401: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1402: @*/
1403: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1404: {
1405:   PetscDualSpaceTransformType trans;
1406:   PetscErrorCode              ierr;

1412:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1413:      This determines their transformation properties. */
1414:   switch (dsp->k)
1415:   {
1416:     case 0: /* H^1 point evaluations */
1417:     trans = IDENTITY_TRANSFORM;break;
1418:     case 1: /* Hcurl preserves tangential edge traces  */
1419:     case 2:
1420:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1421:     case 3: /* Hdiv preserve normal traces */
1422:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1423:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1424:   }
1425:   PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
1426:   return(0);
1427: }

1429: /*@C
1430:   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.

1432:   Input Parameters:
1433: + dsp        - The PetscDualSpace
1434: . fegeom     - The geometry for this cell
1435: . Nq         - The number of function samples
1436: . Nc         - The number of function components
1437: - pointEval  - The function values

1439:   Output Parameter:
1440: . pointEval  - The transformed function values

1442:   Level: advanced

1444:   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.

1446: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1447: @*/
1448: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1449: {
1450:   PetscDualSpaceTransformType trans;
1451:   PetscErrorCode              ierr;

1457:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1458:      This determines their transformation properties. */
1459:   switch (dsp->k)
1460:   {
1461:     case 0: /* H^1 point evaluations */
1462:     trans = IDENTITY_TRANSFORM;break;
1463:     case 1: /* Hcurl preserves tangential edge traces  */
1464:     case 2:
1465:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1466:     case 3: /* Hdiv preserve normal traces */
1467:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1468:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1469:   }
1470:   PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1471:   return(0);
1472: }

1474: /*@C
1475:   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.

1477:   Input Parameters:
1478: + dsp        - The PetscDualSpace
1479: . fegeom     - The geometry for this cell
1480: . Nq         - The number of function gradient samples
1481: . Nc         - The number of function components
1482: - pointEval  - The function gradient values

1484:   Output Parameter:
1485: . pointEval  - The transformed function gradient values

1487:   Level: advanced

1489:   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.

1491: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1492: @*/
1493: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1494: {
1495:   PetscDualSpaceTransformType trans;
1496:   PetscErrorCode              ierr;

1502:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1503:      This determines their transformation properties. */
1504:   switch (dsp->k)
1505:   {
1506:     case 0: /* H^1 point evaluations */
1507:     trans = IDENTITY_TRANSFORM;break;
1508:     case 1: /* Hcurl preserves tangential edge traces  */
1509:     case 2:
1510:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1511:     case 3: /* Hdiv preserve normal traces */
1512:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1513:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1514:   }
1515:   PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1516:   return(0);
1517: }