Actual source code: dualspace.c

petsc-master 2019-10-18
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: beginner

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: intermediate

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:   PetscOptionsBoundedInt("-petscdualspace_degree", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
280:   PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
281:   if (sp->ops->setfromoptions) {
282:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
283:   }
284:   PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
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: intermediate

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: beginner

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: /*@
650:   PetscDualSpaceCreateSection - Create a PetscSection over the reference cell with the layout from this space

652:   Collective on sp

654:   Input Parameters:
655: + sp      - The PetscDualSpace

657:   Output Parameter:
658: . section - The section

660:   Level: advanced

662: .seealso: PetscDualSpaceCreate(), DMPLEX
663: @*/
664: PetscErrorCode PetscDualSpaceCreateSection(PetscDualSpace sp, PetscSection *section)
665: {
666:   DM             dm;
667:   PetscInt       pStart, pEnd, depth, h, offset;
668:   const PetscInt *numDof;

672:   PetscDualSpaceGetDM(sp,&dm);
673:   DMPlexGetChart(dm,&pStart,&pEnd);
674:   PetscSectionCreate(PetscObjectComm((PetscObject)sp),section);
675:   PetscSectionSetChart(*section,pStart,pEnd);
676:   DMPlexGetDepth(dm,&depth);
677:   PetscDualSpaceGetNumDof(sp,&numDof);
678:   for (h = 0; h <= depth; h++) {
679:     PetscInt hStart, hEnd, p, dof;

681:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
682:     dof = numDof[depth - h];
683:     for (p = hStart; p < hEnd; p++) {
684:       PetscSectionSetDof(*section,p,dof);
685:     }
686:   }
687:   PetscSectionSetUp(*section);
688:   for (h = 0, offset = 0; h <= depth; h++) {
689:     PetscInt hStart, hEnd, p, dof;

691:     DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
692:     dof = numDof[depth - h];
693:     for (p = hStart; p < hEnd; p++) {
694:       PetscSectionGetDof(*section,p,&dof);
695:       PetscSectionSetOffset(*section,p,offset);
696:       offset += dof;
697:     }
698:   }
699:   return(0);
700: }

702: /*@
703:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

705:   Collective on sp

707:   Input Parameters:
708: + sp      - The PetscDualSpace
709: . dim     - The spatial dimension
710: - simplex - Flag for simplex, otherwise use a tensor-product cell

712:   Output Parameter:
713: . refdm - The reference cell

715:   Level: intermediate

717: .seealso: PetscDualSpaceCreate(), DMPLEX
718: @*/
719: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
720: {

724:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
725:   return(0);
726: }

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

731:   Input Parameters:
732: + sp      - The PetscDualSpace object
733: . f       - The basis functional index
734: . time    - The time
735: . 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)
736: . numComp - The number of components for the function
737: . func    - The input function
738: - ctx     - A context for the function

740:   Output Parameter:
741: . value   - numComp output values

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

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

748:   Level: beginner

750: .seealso: PetscDualSpaceCreate()
751: @*/
752: 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)
753: {

760:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
761:   return(0);
762: }

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

767:   Input Parameters:
768: + sp        - The PetscDualSpace object
769: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

771:   Output Parameter:
772: . spValue   - The values of all dual space functionals

774:   Level: beginner

776: .seealso: PetscDualSpaceCreate()
777: @*/
778: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
779: {

784:   (*sp->ops->applyall)(sp, pointEval, spValue);
785:   return(0);
786: }

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

791:   Input Parameters:
792: + sp    - The PetscDualSpace object
793: . f     - The basis functional index
794: . time  - The time
795: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
796: . Nc    - The number of components for the function
797: . func  - The input function
798: - ctx   - A context for the function

800:   Output Parameter:
801: . value   - The output value

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

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

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

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

812: where both n and f have Nc components.

814:   Level: beginner

816: .seealso: PetscDualSpaceCreate()
817: @*/
818: 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)
819: {
820:   DM               dm;
821:   PetscQuadrature  n;
822:   const PetscReal *points, *weights;
823:   PetscReal        x[3];
824:   PetscScalar     *val;
825:   PetscInt         dim, dE, qNc, c, Nq, q;
826:   PetscBool        isAffine;
827:   PetscErrorCode   ierr;

832:   PetscDualSpaceGetDM(sp, &dm);
833:   PetscDualSpaceGetFunctional(sp, f, &n);
834:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
835:   if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
836:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
837:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
838:   *value = 0.0;
839:   isAffine = cgeom->isAffine;
840:   dE = cgeom->dimEmbed;
841:   for (q = 0; q < Nq; ++q) {
842:     if (isAffine) {
843:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
844:       (*func)(dE, time, x, Nc, val, ctx);
845:     } else {
846:       (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
847:     }
848:     for (c = 0; c < Nc; ++c) {
849:       *value += val[c]*weights[q*Nc+c];
850:     }
851:   }
852:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
853:   return(0);
854: }

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

859:   Input Parameters:
860: + sp        - The PetscDualSpace object
861: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllPoints()

863:   Output Parameter:
864: . spValue   - The values of all dual space functionals

866:   Level: beginner

868: .seealso: PetscDualSpaceCreate()
869: @*/
870: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
871: {
872:   PetscQuadrature  n;
873:   const PetscReal *points, *weights;
874:   PetscInt         qNc, c, Nq, q, f, spdim, Nc;
875:   PetscInt         offset;
876:   PetscErrorCode   ierr;

882:   PetscDualSpaceGetDimension(sp, &spdim);
883:   PetscDualSpaceGetNumComponents(sp, &Nc);
884:   for (f = 0, offset = 0; f < spdim; f++) {
885:     PetscDualSpaceGetFunctional(sp, f, &n);
886:     PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
887:     if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
888:     spValue[f] = 0.0;
889:     for (q = 0; q < Nq; ++q) {
890:       for (c = 0; c < Nc; ++c) {
891:         spValue[f] += pointEval[offset++]*weights[q*Nc+c];
892:       }
893:     }
894:   }
895:   return(0);
896: }

898: /*@
899:   PetscDualSpaceGetAllPoints - Get all quadrature points from this space

901:   Input Parameter:
902: . sp - The dualspace

904:   Output Parameter:
905: . allPoints - A PetscQuadrature object containing all evaluation points

907:   Level: advanced

909: .seealso: PetscDualSpaceCreate()
910: @*/
911: PetscErrorCode PetscDualSpaceGetAllPoints(PetscDualSpace sp, PetscQuadrature *allPoints)
912: {

918:   if (!sp->allPoints && sp->ops->createallpoints) {
919:     (*sp->ops->createallpoints)(sp,&sp->allPoints);
920:   }
921:   *allPoints = sp->allPoints;
922:   return(0);
923: }

925: /*@
926:   PetscDualSpaceCreateAllPointsDefault - Create all evaluation points by examining functionals

928:   Input Parameter:
929: . sp - The dualspace

931:   Output Parameter:
932: . allPoints - A PetscQuadrature object containing all evaluation points

934:   Level: advanced

936: .seealso: PetscDualSpaceCreate()
937: @*/
938: PetscErrorCode PetscDualSpaceCreateAllPointsDefault(PetscDualSpace sp, PetscQuadrature *allPoints)
939: {
940:   PetscInt        spdim;
941:   PetscInt        numPoints, offset;
942:   PetscReal       *points;
943:   PetscInt        f, dim;
944:   PetscQuadrature q;
945:   PetscErrorCode  ierr;

948:   PetscDualSpaceGetDimension(sp,&spdim);
949:   if (!spdim) {
950:     PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
951:     PetscQuadratureSetData(*allPoints,0,0,0,NULL,NULL);
952:   }
953:   PetscDualSpaceGetFunctional(sp,0,&q);
954:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
955:   for (f = 1; f < spdim; f++) {
956:     PetscInt Np;

958:     PetscDualSpaceGetFunctional(sp,f,&q);
959:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
960:     numPoints += Np;
961:   }
962:   PetscMalloc1(dim*numPoints,&points);
963:   for (f = 0, offset = 0; f < spdim; f++) {
964:     const PetscReal *p;
965:     PetscInt        Np, i;

967:     PetscDualSpaceGetFunctional(sp,f,&q);
968:     PetscQuadratureGetData(q,NULL,NULL,&Np,&p,NULL);
969:     for (i = 0; i < Np * dim; i++) {
970:       points[offset + i] = p[i];
971:     }
972:     offset += Np * dim;
973:   }
974:   PetscQuadratureCreate(PETSC_COMM_SELF,allPoints);
975:   PetscQuadratureSetData(*allPoints,dim,0,numPoints,points,NULL);
976:   return(0);
977: }

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

982:   Input Parameters:
983: + sp    - The PetscDualSpace object
984: . f     - The basis functional index
985: . time  - The time
986: . cgeom - A context with geometric information for this cell, we currently just use the centroid
987: . Nc    - The number of components for the function
988: . func  - The input function
989: - ctx   - A context for the function

991:   Output Parameter:
992: . value - The output value (scalar)

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

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

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

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

1003: where both n and f have Nc components.

1005:   Level: beginner

1007: .seealso: PetscDualSpaceCreate()
1008: @*/
1009: 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)
1010: {
1011:   DM               dm;
1012:   PetscQuadrature  n;
1013:   const PetscReal *points, *weights;
1014:   PetscScalar     *val;
1015:   PetscInt         dimEmbed, qNc, c, Nq, q;
1016:   PetscErrorCode   ierr;

1021:   PetscDualSpaceGetDM(sp, &dm);
1022:   DMGetCoordinateDim(dm, &dimEmbed);
1023:   PetscDualSpaceGetFunctional(sp, f, &n);
1024:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1025:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1026:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1027:   *value = 0.;
1028:   for (q = 0; q < Nq; ++q) {
1029:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1030:     for (c = 0; c < Nc; ++c) {
1031:       *value += val[c]*weights[q*Nc+c];
1032:     }
1033:   }
1034:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1035:   return(0);
1036: }

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

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

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

1048:   Not collective

1050:   Input Parameters:
1051: + sp - the PetscDualSpace object
1052: - height - the height of the mesh point for which the subspace is desired

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

1058:   Level: advanced

1060: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1061: @*/
1062: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1063: {

1069:   *subsp = NULL;
1070:   if (sp->ops->getheightsubspace) {(*sp->ops->getheightsubspace)(sp, height, subsp);}
1071:   return(0);
1072: }

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

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

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

1083:   Not collective

1085:   Input Parameters:
1086: + sp - the PetscDualSpace object
1087: - point - the point (in the dual space's DM) for which the subspace is desired

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

1093:   Level: advanced

1095: .seealso: PetscDualSpace
1096: @*/
1097: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1098: {

1104:   *bdsp = NULL;
1105:   if (sp->ops->getpointsubspace) {
1106:     (*sp->ops->getpointsubspace)(sp,point,bdsp);
1107:   } else if (sp->ops->getheightsubspace) {
1108:     DM       dm;
1109:     DMLabel  label;
1110:     PetscInt dim, depth, height;

1112:     PetscDualSpaceGetDM(sp,&dm);
1113:     DMPlexGetDepth(dm,&dim);
1114:     DMPlexGetDepthLabel(dm,&label);
1115:     DMLabelGetValue(label,point,&depth);
1116:     height = dim - depth;
1117:     (*sp->ops->getheightsubspace)(sp,height,bdsp);
1118:   }
1119:   return(0);
1120: }

1122: /*@C
1123:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1125:   Not collective

1127:   Input Parameter:
1128: . sp - the PetscDualSpace object

1130:   Output Parameters:
1131: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
1132: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation

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

1138:   Level: developer

1140: .seealso: PetscDualSpaceSetSymmetries()
1141: @*/
1142: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1143: {

1150:   if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1151:   return(0);
1152: }

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

1157:   Input Parameter:
1158: . dsp - The PetscDualSpace

1160:   Output Parameter:
1161: . k   - The simplex dimension

1163:   Level: developer

1165:   Note: Currently supported values are
1166: $ 0: These are H_1 methods that only transform coordinates
1167: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1168: $ 2: These are the same as 1
1169: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)

1171: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1172: @*/
1173: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1174: {
1178:   *k = dsp->k;
1179:   return(0);
1180: }

1182: /*@C
1183:   PetscDualSpaceTransform - Transform the function values

1185:   Input Parameters:
1186: + dsp       - The PetscDualSpace
1187: . trans     - The type of transform
1188: . isInverse - Flag to invert the transform
1189: . fegeom    - The cell geometry
1190: . Nv        - The number of function samples
1191: . Nc        - The number of function components
1192: - vals      - The function values

1194:   Output Parameter:
1195: . vals      - The transformed function values

1197:   Level: intermediate

1199: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1200: @*/
1201: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1202: {
1203:   PetscInt dim, v, c;

1209:   dim = dsp->dm->dim;
1210:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1211:   if (Nc == 1 || Nc != dim) return(0);
1212:   switch (trans) {
1213:     case IDENTITY_TRANSFORM: break;
1214:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1215:     if (isInverse) {
1216:       for (v = 0; v < Nv; ++v) {
1217:         switch (dim)
1218:         {
1219:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1220:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1221:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1222:         }
1223:       }
1224:     } else {
1225:       for (v = 0; v < Nv; ++v) {
1226:         switch (dim)
1227:         {
1228:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1229:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1230:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1231:         }
1232:       }
1233:     }
1234:     break;
1235:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1236:     if (isInverse) {
1237:       for (v = 0; v < Nv; ++v) {
1238:         switch (dim)
1239:         {
1240:           case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1241:           case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, 1, &vals[v*Nc], &vals[v*Nc]);break;
1242:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1243:         }
1244:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] *= fegeom->detJ[0];
1245:       }
1246:     } else {
1247:       for (v = 0; v < Nv; ++v) {
1248:         switch (dim)
1249:         {
1250:           case 2: DMPlex_Mult2DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1251:           case 3: DMPlex_Mult3DReal_Internal(fegeom->J, 1, &vals[v*Nc], &vals[v*Nc]);break;
1252:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1253:         }
1254:         for (c = 0; c < Nc; ++c) vals[v*Nc+c] /= fegeom->detJ[0];
1255:       }
1256:     }
1257:     break;
1258:   }
1259:   return(0);
1260: }
1261: /*@C
1262:   PetscDualSpaceTransformGradient - Transform the function gradient values

1264:   Input Parameters:
1265: + dsp       - The PetscDualSpace
1266: . trans     - The type of transform
1267: . isInverse - Flag to invert the transform
1268: . fegeom    - The cell geometry
1269: . Nv        - The number of function gradient samples
1270: . Nc        - The number of function components
1271: - vals      - The function gradient values

1273:   Output Parameter:
1274: . vals      - The transformed function values

1276:   Level: intermediate

1278: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1279: @*/
1280: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1281: {
1282:   PetscInt dim, v, c, d;

1288:   dim = dsp->dm->dim;
1289:   /* Transform gradient */
1290:   for (v = 0; v < Nv; ++v) {
1291:     for (c = 0; c < Nc; ++c) {
1292:       switch (dim)
1293:       {
1294:         case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];
1295:         case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1296:         case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1297:         default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1298:       }
1299:     }
1300:   }
1301:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1302:   if (Nc == 1 || Nc != dim) return(0);
1303:   switch (trans) {
1304:     case IDENTITY_TRANSFORM: break;
1305:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1306:     if (isInverse) {
1307:       for (v = 0; v < Nv; ++v) {
1308:         for (d = 0; d < dim; ++d) {
1309:           switch (dim)
1310:           {
1311:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1312:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1313:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1314:           }
1315:         }
1316:       }
1317:     } else {
1318:       for (v = 0; v < Nv; ++v) {
1319:         for (d = 0; d < dim; ++d) {
1320:           switch (dim)
1321:           {
1322:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1323:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1324:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1325:           }
1326:         }
1327:       }
1328:     }
1329:     break;
1330:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1331:     if (isInverse) {
1332:       for (v = 0; v < Nv; ++v) {
1333:         for (d = 0; d < dim; ++d) {
1334:           switch (dim)
1335:           {
1336:             case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1337:             case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1338:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1339:           }
1340:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1341:         }
1342:       }
1343:     } else {
1344:       for (v = 0; v < Nv; ++v) {
1345:         for (d = 0; d < dim; ++d) {
1346:           switch (dim)
1347:           {
1348:             case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1349:             case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1350:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1351:           }
1352:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1353:         }
1354:       }
1355:     }
1356:     break;
1357:   }
1358:   return(0);
1359: }

1361: /*@C
1362:   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.

1364:   Input Parameters:
1365: + dsp        - The PetscDualSpace
1366: . fegeom     - The geometry for this cell
1367: . Nq         - The number of function samples
1368: . Nc         - The number of function components
1369: - pointEval  - The function values

1371:   Output Parameter:
1372: . pointEval  - The transformed function values

1374:   Level: advanced

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

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

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

1406: /*@C
1407:   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.

1409:   Input Parameters:
1410: + dsp        - The PetscDualSpace
1411: . fegeom     - The geometry for this cell
1412: . Nq         - The number of function samples
1413: . Nc         - The number of function components
1414: - pointEval  - The function values

1416:   Output Parameter:
1417: . pointEval  - The transformed function values

1419:   Level: advanced

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

1423: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1424: @*/
1425: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1426: {
1427:   PetscDualSpaceTransformType trans;
1428:   PetscErrorCode              ierr;

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

1451: /*@C
1452:   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.

1454:   Input Parameters:
1455: + dsp        - The PetscDualSpace
1456: . fegeom     - The geometry for this cell
1457: . Nq         - The number of function gradient samples
1458: . Nc         - The number of function components
1459: - pointEval  - The function gradient values

1461:   Output Parameter:
1462: . pointEval  - The transformed function gradient values

1464:   Level: advanced

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

1468: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
1469: @*/
1470: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
1471: {
1472:   PetscDualSpaceTransformType trans;
1473:   PetscErrorCode              ierr;

1479:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
1480:      This determines their transformation properties. */
1481:   switch (dsp->k)
1482:   {
1483:     case 0: /* H^1 point evaluations */
1484:     trans = IDENTITY_TRANSFORM;break;
1485:     case 1: /* Hcurl preserves tangential edge traces  */
1486:     case 2:
1487:     trans = COVARIANT_PIOLA_TRANSFORM;break;
1488:     case 3: /* Hdiv preserve normal traces */
1489:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
1490:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", dsp->k);
1491:   }
1492:   PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
1493:   return(0);
1494: }