Actual source code: dualspace.c

  1: #include <petsc/private/petscfeimpl.h>
  2: #include <petscdmplex.h>

  4: PetscClassId PETSCDUALSPACE_CLASSID = 0;

  6: PetscLogEvent PETSCDUALSPACE_SetUp;

  8: PetscFunctionList PetscDualSpaceList              = NULL;
  9: PetscBool         PetscDualSpaceRegisterAllCalled = PETSC_FALSE;

 11: const char *const PetscDualSpaceReferenceCells[] = {"SIMPLEX", "TENSOR", "PetscDualSpaceReferenceCell", "PETSCDUALSPACE_REFCELL_", NULL};

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

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

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

 26:   Level: developer

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

 44: /*
 45:   PetscDualSpaceTensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
 46:                                                     Ordering is lexicographic with lowest index as least significant in ordering.
 47:                                                     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}.

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

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

 57:   Level: developer

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

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

 77: /*@C
 78:   PetscDualSpaceRegister - Adds a new PetscDualSpace implementation

 80:   Not Collective

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

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

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

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

104:   Level: advanced

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

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

114:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
115:   return(0);
116: }

118: /*@C
119:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

121:   Collective on sp

123:   Input Parameters:
124: + sp   - The PetscDualSpace object
125: - name - The kind of space

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

130:   Level: intermediate

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

142:   PetscObjectTypeCompare((PetscObject) sp, name, &match);
143:   if (match) return(0);

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

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

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

161:   Not Collective

163:   Input Parameter:
164: . sp  - The PetscDualSpace

166:   Output Parameter:
167: . name - The PetscDualSpace type name

169:   Level: intermediate

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

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

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

194:   PetscDualSpaceGetDimension(sp, &pdim);
195:   PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);
196:   PetscViewerASCIIPushTab(v);
197:   if (sp->k) {
198:     PetscViewerASCIIPrintf(v, "Dual space for %D-forms %swith %D components, size %D\n", PetscAbsInt(sp->k), sp->k < 0 ? "(stored in dual form) ": "", sp->Nc, pdim);
199:   } else {
200:     PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
201:   }
202:   if (sp->ops->view) {(*sp->ops->view)(sp, v);}
203:   PetscViewerGetFormat(v, &format);
204:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
205:     PetscViewerASCIIPushTab(v);
206:     for (f = 0; f < pdim; ++f) {
207:       PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
208:       PetscViewerASCIIPushTab(v);
209:       PetscQuadratureView(sp->functional[f], v);
210:       PetscViewerASCIIPopTab(v);
211:     }
212:     PetscViewerASCIIPopTab(v);
213:   }
214:   PetscViewerASCIIPopTab(v);
215:   return(0);
216: }

218: /*@C
219:    PetscDualSpaceViewFromOptions - View from Options

221:    Collective on PetscDualSpace

223:    Input Parameters:
224: +  A - the PetscDualSpace object
225: .  obj - Optional object, proivides prefix
226: -  name - command line option

228:    Level: intermediate
229: .seealso:  PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate()
230: @*/
231: PetscErrorCode  PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[])
232: {

237:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
238:   return(0);
239: }

241: /*@
242:   PetscDualSpaceView - Views a PetscDualSpace

244:   Collective on sp

246:   Input Parameter:
247: + sp - the PetscDualSpace object to view
248: - v  - the viewer

250:   Level: beginner

252: .seealso PetscDualSpaceDestroy(), PetscDualSpace
253: @*/
254: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
255: {
256:   PetscBool      iascii;

262:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);}
263:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
264:   if (iascii) {PetscDualSpaceView_ASCII(sp, v);}
265:   return(0);
266: }

268: /*@
269:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

271:   Collective on sp

273:   Input Parameter:
274: . sp - the PetscDualSpace object to set options for

276:   Options Database:
277: + -petscdualspace_order <order>      - the approximation order of the space
278: . -petscdualspace_form_degree <deg>  - the form degree, say 0 for point evaluations, or 2 for area integrals
279: . -petscdualspace_components <c>     - the number of components, say d for a vector field
280: . -petscdualspace_refdim <d>         - The spatial dimension of the reference cell
281: - -petscdualspace_refcell <celltype> - Reference cell type name

283:   Level: intermediate

285: .seealso PetscDualSpaceView(), PetscDualSpace, PetscObjectSetFromOptions()
286: @*/
287: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
288: {
289:   PetscDualSpaceReferenceCell refCell = PETSCDUALSPACE_REFCELL_SIMPLEX;
290:   PetscInt                    refDim  = 0;
291:   PetscBool                   flg;
292:   const char                 *defaultType;
293:   char                        name[256];
294:   PetscErrorCode              ierr;

298:   if (!((PetscObject) sp)->type_name) {
299:     defaultType = PETSCDUALSPACELAGRANGE;
300:   } else {
301:     defaultType = ((PetscObject) sp)->type_name;
302:   }
303:   if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}

305:   PetscObjectOptionsBegin((PetscObject) sp);
306:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
307:   if (flg) {
308:     PetscDualSpaceSetType(sp, name);
309:   } else if (!((PetscObject) sp)->type_name) {
310:     PetscDualSpaceSetType(sp, defaultType);
311:   }
312:   PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL,0);
313:   PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);
314:   PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL,1);
315:   if (sp->ops->setfromoptions) {
316:     (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
317:   }
318:   PetscOptionsBoundedInt("-petscdualspace_refdim", "The spatial dimension of the reference cell", "PetscDualSpaceSetReferenceCell", refDim, &refDim, NULL,0);
319:   PetscOptionsEnum("-petscdualspace_refcell", "Reference cell", "PetscDualSpaceSetReferenceCell", PetscDualSpaceReferenceCells, (PetscEnum) refCell, (PetscEnum *) &refCell, &flg);
320:   if (flg) {
321:     DM K;

323:     if (!refDim) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_INCOMP, "Reference cell specified without a dimension. Use -petscdualspace_refdim.");
324:     PetscDualSpaceCreateReferenceCell(sp, refDim, refCell == PETSCDUALSPACE_REFCELL_SIMPLEX ? PETSC_TRUE : PETSC_FALSE, &K);
325:     PetscDualSpaceSetDM(sp, K);
326:     DMDestroy(&K);
327:   }

329:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
330:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
331:   PetscOptionsEnd();
332:   sp->setfromoptionscalled = PETSC_TRUE;
333:   return(0);
334: }

336: /*@
337:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

339:   Collective on sp

341:   Input Parameter:
342: . sp - the PetscDualSpace object to setup

344:   Level: intermediate

346: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace
347: @*/
348: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
349: {

354:   if (sp->setupcalled) return(0);
355:   PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
356:   sp->setupcalled = PETSC_TRUE;
357:   if (sp->ops->setup) {(*sp->ops->setup)(sp);}
358:   PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
359:   if (sp->setfromoptionscalled) {PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");}
360:   return(0);
361: }

363: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
364: {
365:   PetscInt       pStart = -1, pEnd = -1, depth = -1;

369:   if (!dm) return(0);
370:   DMPlexGetChart(dm, &pStart, &pEnd);
371:   DMPlexGetDepth(dm, &depth);

373:   if (sp->pointSpaces) {
374:     PetscInt i;

376:     for (i = 0; i < pEnd - pStart; i++) {
377:       PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
378:     }
379:   }
380:   PetscFree(sp->pointSpaces);

382:   if (sp->heightSpaces) {
383:     PetscInt i;

385:     for (i = 0; i <= depth; i++) {
386:       PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
387:     }
388:   }
389:   PetscFree(sp->heightSpaces);

391:   PetscSectionDestroy(&(sp->pointSection));
392:   PetscQuadratureDestroy(&(sp->intNodes));
393:   VecDestroy(&(sp->intDofValues));
394:   VecDestroy(&(sp->intNodeValues));
395:   MatDestroy(&(sp->intMat));
396:   PetscQuadratureDestroy(&(sp->allNodes));
397:   VecDestroy(&(sp->allDofValues));
398:   VecDestroy(&(sp->allNodeValues));
399:   MatDestroy(&(sp->allMat));
400:   PetscFree(sp->numDof);
401:   return(0);
402: }


405: /*@
406:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

408:   Collective on sp

410:   Input Parameter:
411: . sp - the PetscDualSpace object to destroy

413:   Level: beginner

415: .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate()
416: @*/
417: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
418: {
419:   PetscInt       dim, f;
420:   DM             dm;

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

427:   if (--((PetscObject)(*sp))->refct > 0) {*sp = NULL; return(0);}
428:   ((PetscObject) (*sp))->refct = 0;

430:   PetscDualSpaceGetDimension(*sp, &dim);
431:   dm = (*sp)->dm;

433:   if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
434:   PetscDualSpaceClearDMData_Internal(*sp, dm);

436:   for (f = 0; f < dim; ++f) {
437:     PetscQuadratureDestroy(&(*sp)->functional[f]);
438:   }
439:   PetscFree((*sp)->functional);
440:   DMDestroy(&(*sp)->dm);
441:   PetscHeaderDestroy(sp);
442:   return(0);
443: }

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

448:   Collective

450:   Input Parameter:
451: . comm - The communicator for the PetscDualSpace object

453:   Output Parameter:
454: . sp - The PetscDualSpace object

456:   Level: beginner

458: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
459: @*/
460: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
461: {
462:   PetscDualSpace s;

467:   PetscCitationsRegister(FECitation,&FEcite);
468:   *sp  = NULL;
469:   PetscFEInitializePackage();

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

473:   s->order       = 0;
474:   s->Nc          = 1;
475:   s->k           = 0;
476:   s->spdim       = -1;
477:   s->spintdim    = -1;
478:   s->uniform     = PETSC_TRUE;
479:   s->setupcalled = PETSC_FALSE;

481:   *sp = s;
482:   return(0);
483: }

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

488:   Collective on sp

490:   Input Parameter:
491: . sp - The original PetscDualSpace

493:   Output Parameter:
494: . spNew - The duplicate PetscDualSpace

496:   Level: beginner

498: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
499: @*/
500: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
501: {
502:   DM             dm;
503:   PetscDualSpaceType type;
504:   const char     *name;

510:   PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
511:   PetscObjectGetName((PetscObject) sp,     &name);
512:   PetscObjectSetName((PetscObject) *spNew,  name);
513:   PetscDualSpaceGetType(sp, &type);
514:   PetscDualSpaceSetType(*spNew, type);
515:   PetscDualSpaceGetDM(sp, &dm);
516:   PetscDualSpaceSetDM(*spNew, dm);

518:   (*spNew)->order   = sp->order;
519:   (*spNew)->k       = sp->k;
520:   (*spNew)->Nc      = sp->Nc;
521:   (*spNew)->uniform = sp->uniform;
522:   if (sp->ops->duplicate) {(*sp->ops->duplicate)(sp, *spNew);}
523:   return(0);
524: }

526: /*@
527:   PetscDualSpaceGetDM - Get the DM representing the reference cell

529:   Not collective

531:   Input Parameter:
532: . sp - The PetscDualSpace

534:   Output Parameter:
535: . dm - The reference cell

537:   Level: intermediate

539: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
540: @*/
541: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
542: {
546:   *dm = sp->dm;
547:   return(0);
548: }

550: /*@
551:   PetscDualSpaceSetDM - Get the DM representing the reference cell

553:   Not collective

555:   Input Parameters:
556: + sp - The PetscDualSpace
557: - dm - The reference cell

559:   Level: intermediate

561: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
562: @*/
563: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
564: {

570:   if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change DM after dualspace is set up");
571:   PetscObjectReference((PetscObject) dm);
572:   if (sp->dm && sp->dm != dm) {
573:     PetscDualSpaceClearDMData_Internal(sp, sp->dm);
574:   }
575:   DMDestroy(&sp->dm);
576:   sp->dm = dm;
577:   return(0);
578: }

580: /*@
581:   PetscDualSpaceGetOrder - Get the order of the dual space

583:   Not collective

585:   Input Parameter:
586: . sp - The PetscDualSpace

588:   Output Parameter:
589: . order - The order

591:   Level: intermediate

593: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
594: @*/
595: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
596: {
600:   *order = sp->order;
601:   return(0);
602: }

604: /*@
605:   PetscDualSpaceSetOrder - Set the order of the dual space

607:   Not collective

609:   Input Parameters:
610: + sp - The PetscDualSpace
611: - order - The order

613:   Level: intermediate

615: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
616: @*/
617: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
618: {
621:   if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change order after dualspace is set up");
622:   sp->order = order;
623:   return(0);
624: }

626: /*@
627:   PetscDualSpaceGetNumComponents - Return the number of components for this space

629:   Input Parameter:
630: . sp - The PetscDualSpace

632:   Output Parameter:
633: . Nc - The number of components

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

637:   Level: intermediate

639: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
640: @*/
641: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
642: {
646:   *Nc = sp->Nc;
647:   return(0);
648: }

650: /*@
651:   PetscDualSpaceSetNumComponents - Set the number of components for this space

653:   Input Parameters:
654: + sp - The PetscDualSpace
655: - order - The number of components

657:   Level: intermediate

659: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
660: @*/
661: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
662: {
665:   if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
666:   sp->Nc = Nc;
667:   return(0);
668: }

670: /*@
671:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

673:   Not collective

675:   Input Parameters:
676: + sp - The PetscDualSpace
677: - i  - The basis number

679:   Output Parameter:
680: . functional - The basis functional

682:   Level: intermediate

684: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
685: @*/
686: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
687: {
688:   PetscInt       dim;

694:   PetscDualSpaceGetDimension(sp, &dim);
695:   if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
696:   *functional = sp->functional[i];
697:   return(0);
698: }

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

703:   Not collective

705:   Input Parameter:
706: . sp - The PetscDualSpace

708:   Output Parameter:
709: . dim - The dimension

711:   Level: intermediate

713: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
714: @*/
715: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
716: {

722:   if (sp->spdim < 0) {
723:     PetscSection section;

725:     PetscDualSpaceGetSection(sp, &section);
726:     if (section) {
727:       PetscSectionGetStorageSize(section, &(sp->spdim));
728:     } else sp->spdim = 0;
729:   }
730:   *dim = sp->spdim;
731:   return(0);
732: }

734: /*@
735:   PetscDualSpaceGetInteriorDimension - Get the interior dimension of the dual space, i.e. the number of basis functionals assigned to the interior of the reference domain

737:   Not collective

739:   Input Parameter:
740: . sp - The PetscDualSpace

742:   Output Parameter:
743: . dim - The dimension

745:   Level: intermediate

747: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
748: @*/
749: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
750: {

756:   if (sp->spintdim < 0) {
757:     PetscSection section;

759:     PetscDualSpaceGetSection(sp, &section);
760:     if (section) {
761:       PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
762:     } else sp->spintdim = 0;
763:   }
764:   *intdim = sp->spintdim;
765:   return(0);
766: }

768: /*@
769:    PetscDualSpaceGetUniform - Whether this dual space is uniform

771:    Not collective

773:    Input Parameters:
774: .  sp - A dual space

776:    Output Parameters:
777: .  uniform - PETSC_TRUE if (a) the dual space is the same for each point in a stratum of the reference DMPlex, and
778:              (b) every symmetry of each point in the reference DMPlex is also a symmetry of the point's dual space.


781:    Level: advanced

783:    Note: all of the usual spaces on simplex or tensor-product elements will be uniform, only reference cells
784:    with non-uniform strata (like trianguar-prisms) or anisotropic hp dual spaces will not be uniform.

786: .seealso: PetscDualSpaceGetPointSubspace(), PetscDualSpaceGetSymmetries()
787: @*/
788: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
789: {
793:   *uniform = sp->uniform;
794:   return(0);
795: }


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

801:   Not collective

803:   Input Parameter:
804: . sp - The PetscDualSpace

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

809:   Level: intermediate

811: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
812: @*/
813: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
814: {

820:   if (!sp->uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform space does not have a fixed number of dofs for each height");
821:   if (!sp->numDof) {
822:     DM       dm;
823:     PetscInt depth, d;
824:     PetscSection section;

826:     PetscDualSpaceGetDM(sp, &dm);
827:     DMPlexGetDepth(dm, &depth);
828:     PetscCalloc1(depth+1,&(sp->numDof));
829:     PetscDualSpaceGetSection(sp, &section);
830:     for (d = 0; d <= depth; d++) {
831:       PetscInt dStart, dEnd;

833:       DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
834:       if (dEnd <= dStart) continue;
835:       PetscSectionGetDof(section, dStart, &(sp->numDof[d]));

837:     }
838:   }
839:   *numDof = sp->numDof;
840:   if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
841:   return(0);
842: }

844: /* create the section of the right size and set a permutation for topological ordering */
845: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
846: {
847:   DM             dm;
848:   PetscInt       pStart, pEnd, cStart, cEnd, c, depth, count, i;
849:   PetscInt       *seen, *perm;
850:   PetscSection   section;

854:   dm = sp->dm;
855:   PetscSectionCreate(PETSC_COMM_SELF, &section);
856:   DMPlexGetChart(dm, &pStart, &pEnd);
857:   PetscSectionSetChart(section, pStart, pEnd);
858:   PetscCalloc1(pEnd - pStart, &seen);
859:   PetscMalloc1(pEnd - pStart, &perm);
860:   DMPlexGetDepth(dm, &depth);
861:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
862:   for (c = cStart, count = 0; c < cEnd; c++) {
863:     PetscInt closureSize = -1, e;
864:     PetscInt *closure = NULL;

866:     perm[count++] = c;
867:     seen[c-pStart] = 1;
868:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
869:     for (e = 0; e < closureSize; e++) {
870:       PetscInt point = closure[2*e];

872:       if (seen[point-pStart]) continue;
873:       perm[count++] = point;
874:       seen[point-pStart] = 1;
875:     }
876:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
877:   }
878:   if (count != pEnd - pStart) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad topological ordering");
879:   for (i = 0; i < pEnd - pStart; i++) if (perm[i] != i) break;
880:   if (i < pEnd - pStart) {
881:     IS permIS;

883:     ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
884:     ISSetPermutation(permIS);
885:     PetscSectionSetPermutation(section, permIS);
886:     ISDestroy(&permIS);
887:   } else {
888:     PetscFree(perm);
889:   }
890:   PetscFree(seen);
891:   *topSection = section;
892:   return(0);
893: }

895: /* mark boundary points and set up */
896: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
897: {
898:   DM             dm;
899:   DMLabel        boundary;
900:   PetscInt       pStart, pEnd, p;

904:   dm = sp->dm;
905:   DMLabelCreate(PETSC_COMM_SELF,"boundary",&boundary);
906:   PetscDualSpaceGetDM(sp,&dm);
907:   DMPlexMarkBoundaryFaces(dm,1,boundary);
908:   DMPlexLabelComplete(dm,boundary);
909:   DMPlexGetChart(dm, &pStart, &pEnd);
910:   for (p = pStart; p < pEnd; p++) {
911:     PetscInt bval;

913:     DMLabelGetValue(boundary, p, &bval);
914:     if (bval == 1) {
915:       PetscInt dof;

917:       PetscSectionGetDof(section, p, &dof);
918:       PetscSectionSetConstraintDof(section, p, dof);
919:     }
920:   }
921:   DMLabelDestroy(&boundary);
922:   PetscSectionSetUp(section);
923:   return(0);
924: }

926: /*@
927:   PetscDualSpaceGetSection - Create a PetscSection over the reference cell with the layout from this space

929:   Collective on sp

931:   Input Parameters:
932: . sp      - The PetscDualSpace

934:   Output Parameter:
935: . section - The section

937:   Level: advanced

939: .seealso: PetscDualSpaceCreate(), DMPLEX
940: @*/
941: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
942: {
943:   PetscInt       pStart, pEnd, p;

947:   if (!sp->pointSection) {
948:     /* mark the boundary */
949:     PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
950:     DMPlexGetChart(sp->dm,&pStart,&pEnd);
951:     for (p = pStart; p < pEnd; p++) {
952:       PetscDualSpace psp;

954:       PetscDualSpaceGetPointSubspace(sp, p, &psp);
955:       if (psp) {
956:         PetscInt dof;

958:         PetscDualSpaceGetInteriorDimension(psp, &dof);
959:         PetscSectionSetDof(sp->pointSection,p,dof);
960:       }
961:     }
962:     PetscDualSpaceSectionSetUp_Internal(sp,sp->pointSection);
963:   }
964:   *section = sp->pointSection;
965:   return(0);
966: }

968: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
969:  * have one cell */
970: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
971: {
972:   PetscReal *sv0, *v0, *J;
973:   PetscSection section;
974:   PetscInt     dim, s, k;
975:   DM             dm;

979:   PetscDualSpaceGetDM(sp, &dm);
980:   DMGetDimension(dm, &dim);
981:   PetscDualSpaceGetSection(sp, &section);
982:   PetscMalloc3(dim, &v0, dim, &sv0, dim*dim, &J);
983:   PetscDualSpaceGetFormDegree(sp, &k);
984:   for (s = sStart; s < sEnd; s++) {
985:     PetscReal detJ, hdetJ;
986:     PetscDualSpace ssp;
987:     PetscInt dof, off, f, sdim;
988:     PetscInt i, j;
989:     DM sdm;

991:     PetscDualSpaceGetPointSubspace(sp, s, &ssp);
992:     if (!ssp) continue;
993:     PetscSectionGetDof(section, s, &dof);
994:     PetscSectionGetOffset(section, s, &off);
995:     /* get the first vertex of the reference cell */
996:     PetscDualSpaceGetDM(ssp, &sdm);
997:     DMGetDimension(sdm, &sdim);
998:     DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);
999:     DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);
1000:     /* compactify Jacobian */
1001:     for (i = 0; i < dim; i++) for (j = 0; j < sdim; j++) J[i* sdim + j] = J[i * dim + j];
1002:     for (f = 0; f < dof; f++) {
1003:       PetscQuadrature fn;

1005:       PetscDualSpaceGetFunctional(ssp, f, &fn);
1006:       PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off+f]));
1007:     }
1008:   }
1009:   PetscFree3(v0, sv0, J);
1010:   return(0);
1011: }

1013: /*@
1014:   PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell

1016:   Collective on sp

1018:   Input Parameters:
1019: + sp      - The PetscDualSpace
1020: . dim     - The spatial dimension
1021: - simplex - Flag for simplex, otherwise use a tensor-product cell

1023:   Output Parameter:
1024: . refdm - The reference cell

1026:   Level: intermediate

1028: .seealso: PetscDualSpaceCreate(), DMPLEX
1029: @*/
1030: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1031: {

1035:   DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1036:   return(0);
1037: }

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

1042:   Input Parameters:
1043: + sp      - The PetscDualSpace object
1044: . f       - The basis functional index
1045: . time    - The time
1046: . 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)
1047: . numComp - The number of components for the function
1048: . func    - The input function
1049: - ctx     - A context for the function

1051:   Output Parameter:
1052: . value   - numComp output values

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

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

1059:   Level: beginner

1061: .seealso: PetscDualSpaceCreate()
1062: @*/
1063: 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)
1064: {

1071:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
1072:   return(0);
1073: }

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

1078:   Input Parameters:
1079: + sp        - The PetscDualSpace object
1080: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()

1082:   Output Parameter:
1083: . spValue   - The values of all dual space functionals

1085:   Level: beginner

1087: .seealso: PetscDualSpaceCreate()
1088: @*/
1089: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1090: {

1095:   (*sp->ops->applyall)(sp, pointEval, spValue);
1096:   return(0);
1097: }

1099: /*@C
1100:   PetscDualSpaceApplyInterior - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()

1102:   Input Parameters:
1103: + sp        - The PetscDualSpace object
1104: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()

1106:   Output Parameter:
1107: . spValue   - The values of interior dual space functionals

1109:   Level: beginner

1111: .seealso: PetscDualSpaceCreate()
1112: @*/
1113: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1114: {

1119:   (*sp->ops->applyint)(sp, pointEval, spValue);
1120:   return(0);
1121: }

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

1126:   Input Parameters:
1127: + sp    - The PetscDualSpace object
1128: . f     - The basis functional index
1129: . time  - The time
1130: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1131: . Nc    - The number of components for the function
1132: . func  - The input function
1133: - ctx   - A context for the function

1135:   Output Parameter:
1136: . value   - The output value

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

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

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

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

1147: where both n and f have Nc components.

1149:   Level: beginner

1151: .seealso: PetscDualSpaceCreate()
1152: @*/
1153: 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)
1154: {
1155:   DM               dm;
1156:   PetscQuadrature  n;
1157:   const PetscReal *points, *weights;
1158:   PetscReal        x[3];
1159:   PetscScalar     *val;
1160:   PetscInt         dim, dE, qNc, c, Nq, q;
1161:   PetscBool        isAffine;
1162:   PetscErrorCode   ierr;

1167:   PetscDualSpaceGetDM(sp, &dm);
1168:   PetscDualSpaceGetFunctional(sp, f, &n);
1169:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
1170:   if (dim != cgeom->dim) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature spatial dimension %D != cell geometry dimension %D", dim, cgeom->dim);
1171:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1172:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1173:   *value = 0.0;
1174:   isAffine = cgeom->isAffine;
1175:   dE = cgeom->dimEmbed;
1176:   for (q = 0; q < Nq; ++q) {
1177:     if (isAffine) {
1178:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q*dim], x);
1179:       (*func)(dE, time, x, Nc, val, ctx);
1180:     } else {
1181:       (*func)(dE, time, &cgeom->v[dE*q], Nc, val, ctx);
1182:     }
1183:     for (c = 0; c < Nc; ++c) {
1184:       *value += val[c]*weights[q*Nc+c];
1185:     }
1186:   }
1187:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1188:   return(0);
1189: }

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

1194:   Input Parameters:
1195: + sp        - The PetscDualSpace object
1196: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()

1198:   Output Parameter:
1199: . spValue   - The values of all dual space functionals

1201:   Level: beginner

1203: .seealso: PetscDualSpaceCreate()
1204: @*/
1205: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1206: {
1207:   Vec              pointValues, dofValues;
1208:   Mat              allMat;
1209:   PetscErrorCode   ierr;

1215:   PetscDualSpaceGetAllData(sp, NULL, &allMat);
1216:   if (!(sp->allNodeValues)) {
1217:     MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1218:   }
1219:   pointValues = sp->allNodeValues;
1220:   if (!(sp->allDofValues)) {
1221:     MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1222:   }
1223:   dofValues = sp->allDofValues;
1224:   VecPlaceArray(pointValues, pointEval);
1225:   VecPlaceArray(dofValues, spValue);
1226:   MatMult(allMat, pointValues, dofValues);
1227:   VecResetArray(dofValues);
1228:   VecResetArray(pointValues);
1229:   return(0);
1230: }

1232: /*@C
1233:   PetscDualSpaceApplyInteriorDefault - Apply interior functionals from the dual space basis to the result of an evaluation at the points returned by PetscDualSpaceGetInteriorData()

1235:   Input Parameters:
1236: + sp        - The PetscDualSpace object
1237: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()

1239:   Output Parameter:
1240: . spValue   - The values of interior dual space functionals

1242:   Level: beginner

1244: .seealso: PetscDualSpaceCreate()
1245: @*/
1246: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1247: {
1248:   Vec              pointValues, dofValues;
1249:   Mat              intMat;
1250:   PetscErrorCode   ierr;

1256:   PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1257:   if (!(sp->intNodeValues)) {
1258:     MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1259:   }
1260:   pointValues = sp->intNodeValues;
1261:   if (!(sp->intDofValues)) {
1262:     MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1263:   }
1264:   dofValues = sp->intDofValues;
1265:   VecPlaceArray(pointValues, pointEval);
1266:   VecPlaceArray(dofValues, spValue);
1267:   MatMult(intMat, pointValues, dofValues);
1268:   VecResetArray(dofValues);
1269:   VecResetArray(pointValues);
1270:   return(0);
1271: }

1273: /*@
1274:   PetscDualSpaceGetAllData - Get all quadrature nodes from this space, and the matrix that sends quadrature node values to degree-of-freedom values

1276:   Input Parameter:
1277: . sp - The dualspace

1279:   Output Parameter:
1280: + allNodes - A PetscQuadrature object containing all evaluation nodes
1281: - allMat - A Mat for the node-to-dof transformation

1283:   Level: advanced

1285: .seealso: PetscDualSpaceCreate()
1286: @*/
1287: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1288: {

1295:   if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1296:     PetscQuadrature qpoints;
1297:     Mat amat;

1299:     (*sp->ops->createalldata)(sp,&qpoints,&amat);
1300:     PetscQuadratureDestroy(&(sp->allNodes));
1301:     MatDestroy(&(sp->allMat));
1302:     sp->allNodes = qpoints;
1303:     sp->allMat = amat;
1304:   }
1305:   if (allNodes) *allNodes = sp->allNodes;
1306:   if (allMat) *allMat = sp->allMat;
1307:   return(0);
1308: }

1310: /*@
1311:   PetscDualSpaceCreateAllDataDefault - Create all evaluation nodes and the node-to-dof matrix by examining functionals

1313:   Input Parameter:
1314: . sp - The dualspace

1316:   Output Parameter:
1317: + allNodes - A PetscQuadrature object containing all evaluation nodes
1318: - allMat - A Mat for the node-to-dof transformation

1320:   Level: advanced

1322: .seealso: PetscDualSpaceCreate()
1323: @*/
1324: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1325: {
1326:   PetscInt        spdim;
1327:   PetscInt        numPoints, offset;
1328:   PetscReal       *points;
1329:   PetscInt        f, dim;
1330:   PetscInt        Nc, nrows, ncols;
1331:   PetscInt        maxNumPoints;
1332:   PetscQuadrature q;
1333:   Mat             A;
1334:   PetscErrorCode  ierr;

1337:   PetscDualSpaceGetNumComponents(sp, &Nc);
1338:   PetscDualSpaceGetDimension(sp,&spdim);
1339:   if (!spdim) {
1340:     PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1341:     PetscQuadratureSetData(*allNodes,0,0,0,NULL,NULL);
1342:   }
1343:   nrows = spdim;
1344:   PetscDualSpaceGetFunctional(sp,0,&q);
1345:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
1346:   maxNumPoints = numPoints;
1347:   for (f = 1; f < spdim; f++) {
1348:     PetscInt Np;

1350:     PetscDualSpaceGetFunctional(sp,f,&q);
1351:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1352:     numPoints += Np;
1353:     maxNumPoints = PetscMax(maxNumPoints,Np);
1354:   }
1355:   ncols = numPoints * Nc;
1356:   PetscMalloc1(dim*numPoints,&points);
1357:   MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1358:   for (f = 0, offset = 0; f < spdim; f++) {
1359:     const PetscReal *p, *w;
1360:     PetscInt        Np, i;
1361:     PetscInt        fnc;

1363:     PetscDualSpaceGetFunctional(sp,f,&q);
1364:     PetscQuadratureGetData(q,NULL,&fnc,&Np,&p,&w);
1365:     if (fnc != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "functional component mismatch");
1366:     for (i = 0; i < Np * dim; i++) {
1367:       points[offset* dim + i] = p[i];
1368:     }
1369:     for (i = 0; i < Np * Nc; i++) {
1370:       MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1371:     }
1372:     offset += Np;
1373:   }
1374:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1375:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1376:   PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1377:   PetscQuadratureSetData(*allNodes,dim,0,numPoints,points,NULL);
1378:   *allMat = A;
1379:   return(0);
1380: }

1382: /*@
1383:   PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1384:   this space, as well as the matrix that computes the degrees of freedom from the quadrature values.  Degrees of
1385:   freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the
1386:   reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by
1387:   PetscDualSpaceGetSection()).

1389:   Input Parameter:
1390: . sp - The dualspace

1392:   Output Parameter:
1393: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1394: - intMat   - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1395:              the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1396:              npoints is the number of points in intNodes and nc is PetscDualSpaceGetNumComponents().

1398:   Level: advanced

1400: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetDimension(), PetscDualSpaceGetNumComponents(), PetscQuadratureGetData()
1401: @*/
1402: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1403: {

1410:   if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1411:     PetscQuadrature qpoints;
1412:     Mat imat;

1414:     (*sp->ops->createintdata)(sp,&qpoints,&imat);
1415:     PetscQuadratureDestroy(&(sp->intNodes));
1416:     MatDestroy(&(sp->intMat));
1417:     sp->intNodes = qpoints;
1418:     sp->intMat = imat;
1419:   }
1420:   if (intNodes) *intNodes = sp->intNodes;
1421:   if (intMat) *intMat = sp->intMat;
1422:   return(0);
1423: }

1425: /*@
1426:   PetscDualSpaceCreateInteriorDataDefault - Create quadrature points by examining interior functionals and create the matrix mapping quadrature point values to interior dual space values

1428:   Input Parameter:
1429: . sp - The dualspace

1431:   Output Parameter:
1432: + intNodes - A PetscQuadrature object containing all evaluation points needed to evaluate interior degrees of freedom
1433: - intMat    - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1434:               the size of the constrained layout (PetscSectionGetConstrainStorageSize()) of the dual space section,
1435:               npoints is the number of points in allNodes and nc is PetscDualSpaceGetNumComponents().

1437:   Level: advanced

1439: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetInteriorData()
1440: @*/
1441: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1442: {
1443:   DM              dm;
1444:   PetscInt        spdim0;
1445:   PetscInt        Nc;
1446:   PetscInt        pStart, pEnd, p, f;
1447:   PetscSection    section;
1448:   PetscInt        numPoints, offset, matoffset;
1449:   PetscReal       *points;
1450:   PetscInt        dim;
1451:   PetscInt        *nnz;
1452:   PetscQuadrature q;
1453:   Mat             imat;
1454:   PetscErrorCode  ierr;

1458:   PetscDualSpaceGetSection(sp, &section);
1459:   PetscSectionGetConstrainedStorageSize(section, &spdim0);
1460:   if (!spdim0) {
1461:     *intNodes = NULL;
1462:     *intMat = NULL;
1463:     return(0);
1464:   }
1465:   PetscDualSpaceGetNumComponents(sp, &Nc);
1466:   PetscSectionGetChart(section, &pStart, &pEnd);
1467:   PetscDualSpaceGetDM(sp, &dm);
1468:   DMGetDimension(dm, &dim);
1469:   PetscMalloc1(spdim0, &nnz);
1470:   for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1471:     PetscInt dof, cdof, off, d;

1473:     PetscSectionGetDof(section, p, &dof);
1474:     PetscSectionGetConstraintDof(section, p, &cdof);
1475:     if (!(dof - cdof)) continue;
1476:     PetscSectionGetOffset(section, p, &off);
1477:     for (d = 0; d < dof; d++, off++, f++) {
1478:       PetscInt Np;

1480:       PetscDualSpaceGetFunctional(sp,off,&q);
1481:       PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1482:       nnz[f] = Np * Nc;
1483:       numPoints += Np;
1484:     }
1485:   }
1486:   MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1487:   PetscFree(nnz);
1488:   PetscMalloc1(dim*numPoints,&points);
1489:   for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1490:     PetscInt dof, cdof, off, d;

1492:     PetscSectionGetDof(section, p, &dof);
1493:     PetscSectionGetConstraintDof(section, p, &cdof);
1494:     if (!(dof - cdof)) continue;
1495:     PetscSectionGetOffset(section, p, &off);
1496:     for (d = 0; d < dof; d++, off++, f++) {
1497:       const PetscReal *p;
1498:       const PetscReal *w;
1499:       PetscInt        Np, i;

1501:       PetscDualSpaceGetFunctional(sp,off,&q);
1502:       PetscQuadratureGetData(q,NULL,NULL,&Np,&p,&w);
1503:       for (i = 0; i < Np * dim; i++) {
1504:         points[offset + i] = p[i];
1505:       }
1506:       for (i = 0; i < Np * Nc; i++) {
1507:         MatSetValue(imat, f, matoffset + i, w[i],INSERT_VALUES);
1508:       }
1509:       offset += Np * dim;
1510:       matoffset += Np * Nc;
1511:     }
1512:   }
1513:   PetscQuadratureCreate(PETSC_COMM_SELF,intNodes);
1514:   PetscQuadratureSetData(*intNodes,dim,0,numPoints,points,NULL);
1515:   MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1516:   MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1517:   *intMat = imat;
1518:   return(0);
1519: }

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

1524:   Input Parameters:
1525: + sp    - The PetscDualSpace object
1526: . f     - The basis functional index
1527: . time  - The time
1528: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1529: . Nc    - The number of components for the function
1530: . func  - The input function
1531: - ctx   - A context for the function

1533:   Output Parameter:
1534: . value - The output value (scalar)

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

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

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

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

1545: where both n and f have Nc components.

1547:   Level: beginner

1549: .seealso: PetscDualSpaceCreate()
1550: @*/
1551: 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)
1552: {
1553:   DM               dm;
1554:   PetscQuadrature  n;
1555:   const PetscReal *points, *weights;
1556:   PetscScalar     *val;
1557:   PetscInt         dimEmbed, qNc, c, Nq, q;
1558:   PetscErrorCode   ierr;

1563:   PetscDualSpaceGetDM(sp, &dm);
1564:   DMGetCoordinateDim(dm, &dimEmbed);
1565:   PetscDualSpaceGetFunctional(sp, f, &n);
1566:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1567:   if (qNc != Nc) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_SIZ, "The quadrature components %D != function components %D", qNc, Nc);
1568:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1569:   *value = 0.;
1570:   for (q = 0; q < Nq; ++q) {
1571:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1572:     for (c = 0; c < Nc; ++c) {
1573:       *value += val[c]*weights[q*Nc+c];
1574:     }
1575:   }
1576:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1577:   return(0);
1578: }

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

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

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

1590:   Not collective

1592:   Input Parameters:
1593: + sp - the PetscDualSpace object
1594: - height - the height of the mesh point for which the subspace is desired

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

1600:   Level: advanced

1602: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1603: @*/
1604: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1605: {
1606:   PetscInt       depth = -1, cStart, cEnd;
1607:   DM             dm;

1613:   if (!(sp->uniform)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "A non-uniform dual space does not have a single dual space at each height");
1614:   *subsp = NULL;
1615:   dm = sp->dm;
1616:   DMPlexGetDepth(dm, &depth);
1617:   if (height < 0 || height > depth) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid height");
1618:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1619:   if (height == 0 && cEnd == cStart + 1) {
1620:     *subsp = sp;
1621:     return(0);
1622:   }
1623:   if (!sp->heightSpaces) {
1624:     PetscInt h;
1625:     PetscCalloc1(depth+1, &(sp->heightSpaces));

1627:     for (h = 0; h <= depth; h++) {
1628:       if (h == 0 && cEnd == cStart + 1) continue;
1629:       if (sp->ops->createheightsubspace) {(*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));}
1630:       else if (sp->pointSpaces) {
1631:         PetscInt hStart, hEnd;

1633:         DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
1634:         if (hEnd > hStart) {
1635:           const char *name;

1637:           PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1638:           if (sp->pointSpaces[hStart]) {
1639:             PetscObjectGetName((PetscObject) sp,                     &name);
1640:             PetscObjectSetName((PetscObject) sp->pointSpaces[hStart], name);
1641:           }
1642:           sp->heightSpaces[h] = sp->pointSpaces[hStart];
1643:         }
1644:       }
1645:     }
1646:   }
1647:   *subsp = sp->heightSpaces[height];
1648:   return(0);
1649: }

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

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

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

1660:   Not collective

1662:   Input Parameters:
1663: + sp - the PetscDualSpace object
1664: - point - the point (in the dual space's DM) for which the subspace is desired

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

1670:   Level: advanced

1672: .seealso: PetscDualSpace
1673: @*/
1674: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1675: {
1676:   PetscInt       pStart = 0, pEnd = 0, cStart, cEnd;
1677:   DM             dm;

1683:   *bdsp = NULL;
1684:   dm = sp->dm;
1685:   DMPlexGetChart(dm, &pStart, &pEnd);
1686:   if (point < pStart || point > pEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point");
1687:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1688:   if (point == cStart && cEnd == cStart + 1) { /* the dual space is only equivalent to the dual space on a cell if the reference mesh has just one cell */
1689:     *bdsp = sp;
1690:     return(0);
1691:   }
1692:   if (!sp->pointSpaces) {
1693:     PetscInt p;
1694:     PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));

1696:     for (p = 0; p < pEnd - pStart; p++) {
1697:       if (p + pStart == cStart && cEnd == cStart + 1) continue;
1698:       if (sp->ops->createpointsubspace) {(*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));}
1699:       else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1700:         PetscInt dim, depth, height;
1701:         DMLabel  label;

1703:         DMPlexGetDepth(dm,&dim);
1704:         DMPlexGetDepthLabel(dm,&label);
1705:         DMLabelGetValue(label,p+pStart,&depth);
1706:         height = dim - depth;
1707:         PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1708:         PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1709:       }
1710:     }
1711:   }
1712:   *bdsp = sp->pointSpaces[point - pStart];
1713:   return(0);
1714: }

1716: /*@C
1717:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1719:   Not collective

1721:   Input Parameter:
1722: . sp - the PetscDualSpace object

1724:   Output Parameters:
1725: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1726: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation

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

1732:   Level: developer

1734: @*/
1735: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1736: {

1743:   if (sp->ops->getsymmetries) {(sp->ops->getsymmetries)(sp,perms,flips);}
1744:   return(0);
1745: }

1747: /*@
1748:   PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1749:   dual space's functionals.

1751:   Input Parameter:
1752: . dsp - The PetscDualSpace

1754:   Output Parameter:
1755: . k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1756:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1757:         the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz).
1758:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1759:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1760:         but are stored as 1-forms.

1762:   Level: developer

1764: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1765: @*/
1766: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1767: {
1771:   *k = dsp->k;
1772:   return(0);
1773: }

1775: /*@
1776:   PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1777:   dual space's functionals.

1779:   Input Parameter:
1780: + dsp - The PetscDualSpace
1781: - k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1782:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1783:         the 1-form basis in 3-D is (dx, dy, dz), and the 2-form basis in 3-D is (dx wedge dy, dx wedge dz, dy wedge dz).
1784:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1785:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1786:         but are stored as 1-forms.

1788:   Level: developer

1790: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1791: @*/
1792: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1793: {
1794:   PetscInt dim;

1798:   if (dsp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of components after dualspace is set up");
1799:   dim = dsp->dm->dim;
1800:   if (k < -dim || k > dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported %D-form on %D-dimensional reference cell", PetscAbsInt(k), dim);
1801:   dsp->k = k;
1802:   return(0);
1803: }

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

1808:   Input Parameter:
1809: . dsp - The PetscDualSpace

1811:   Output Parameter:
1812: . k   - The simplex dimension

1814:   Level: developer

1816:   Note: Currently supported values are
1817: $ 0: These are H_1 methods that only transform coordinates
1818: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1819: $ 2: These are the same as 1
1820: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)

1822: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1823: @*/
1824: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1825: {
1826:   PetscInt dim;

1831:   dim = dsp->dm->dim;
1832:   if (!dsp->k) *k = IDENTITY_TRANSFORM;
1833:   else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1834:   else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1835:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1836:   return(0);
1837: }

1839: /*@C
1840:   PetscDualSpaceTransform - Transform the function values

1842:   Input Parameters:
1843: + dsp       - The PetscDualSpace
1844: . trans     - The type of transform
1845: . isInverse - Flag to invert the transform
1846: . fegeom    - The cell geometry
1847: . Nv        - The number of function samples
1848: . Nc        - The number of function components
1849: - vals      - The function values

1851:   Output Parameter:
1852: . vals      - The transformed function values

1854:   Level: intermediate

1856:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1858: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpaceTransformHessian(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1859: @*/
1860: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1861: {
1862:   PetscReal Jstar[9] = {0};
1863:   PetscInt dim, v, c, Nk;

1870:   /* TODO: not handling dimEmbed != dim right now */
1871:   dim = dsp->dm->dim;
1872:   /* No change needed for 0-forms */
1873:   if (!dsp->k) return(0);
1874:   PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1875:   /* TODO: use fegeom->isAffine */
1876:   PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1877:   for (v = 0; v < Nv; ++v) {
1878:     switch (Nk) {
1879:     case 1:
1880:       for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0];
1881:       break;
1882:     case 2:
1883:       for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1884:       break;
1885:     case 3:
1886:       for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1887:       break;
1888:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk);
1889:     }
1890:   }
1891:   return(0);
1892: }

1894: /*@C
1895:   PetscDualSpaceTransformGradient - Transform the function gradient values

1897:   Input Parameters:
1898: + dsp       - The PetscDualSpace
1899: . trans     - The type of transform
1900: . isInverse - Flag to invert the transform
1901: . fegeom    - The cell geometry
1902: . Nv        - The number of function gradient samples
1903: . Nc        - The number of function components
1904: - vals      - The function gradient values

1906:   Output Parameter:
1907: . vals      - The transformed function gradient values

1909:   Level: intermediate

1911:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

1913: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1914: @*/
1915: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1916: {
1917:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1918:   PetscInt       v, c, d;

1924: #ifdef PETSC_USE_DEBUG
1925:   if (dE <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE);
1926: #endif
1927:   /* Transform gradient */
1928:   if (dim == dE) {
1929:     for (v = 0; v < Nv; ++v) {
1930:       for (c = 0; c < Nc; ++c) {
1931:         switch (dim)
1932:         {
1933:           case 1: vals[(v*Nc+c)*dim] *= fegeom->invJ[0];break;
1934:           case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1935:           case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v*Nc+c)*dim], &vals[(v*Nc+c)*dim]);break;
1936:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1937:         }
1938:       }
1939:     }
1940:   } else {
1941:     for (v = 0; v < Nv; ++v) {
1942:       for (c = 0; c < Nc; ++c) {
1943:         DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v*Nc+c)*dE], &vals[(v*Nc+c)*dE]);
1944:       }
1945:     }
1946:   }
1947:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1948:   if (Nc == 1 || Nc != dim) return(0);
1949:   switch (trans) {
1950:     case IDENTITY_TRANSFORM: break;
1951:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1952:     if (isInverse) {
1953:       for (v = 0; v < Nv; ++v) {
1954:         for (d = 0; d < dim; ++d) {
1955:           switch (dim)
1956:           {
1957:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1958:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1959:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1960:           }
1961:         }
1962:       }
1963:     } else {
1964:       for (v = 0; v < Nv; ++v) {
1965:         for (d = 0; d < dim; ++d) {
1966:           switch (dim)
1967:           {
1968:             case 2: DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1969:             case 3: DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1970:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1971:           }
1972:         }
1973:       }
1974:     }
1975:     break;
1976:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1977:     if (isInverse) {
1978:       for (v = 0; v < Nv; ++v) {
1979:         for (d = 0; d < dim; ++d) {
1980:           switch (dim)
1981:           {
1982:             case 2: DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1983:             case 3: DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1984:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1985:           }
1986:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] *= fegeom->detJ[0];
1987:         }
1988:       }
1989:     } else {
1990:       for (v = 0; v < Nv; ++v) {
1991:         for (d = 0; d < dim; ++d) {
1992:           switch (dim)
1993:           {
1994:             case 2: DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1995:             case 3: DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v*Nc*dim+d], &vals[v*Nc*dim+d]);break;
1996:             default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
1997:           }
1998:           for (c = 0; c < Nc; ++c) vals[(v*Nc+c)*dim+d] /= fegeom->detJ[0];
1999:         }
2000:       }
2001:     }
2002:     break;
2003:   }
2004:   return(0);
2005: }

2007: /*@C
2008:   PetscDualSpaceTransformHessian - Transform the function Hessian values

2010:   Input Parameters:
2011: + dsp       - The PetscDualSpace
2012: . trans     - The type of transform
2013: . isInverse - Flag to invert the transform
2014: . fegeom    - The cell geometry
2015: . Nv        - The number of function Hessian samples
2016: . Nc        - The number of function components
2017: - vals      - The function gradient values

2019:   Output Parameter:
2020: . vals      - The transformed function Hessian values

2022:   Level: intermediate

2024:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

2026: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
2027: @*/
2028: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
2029: {
2030:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
2031:   PetscInt       v, c;

2037: #ifdef PETSC_USE_DEBUG
2038:   if (dE <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid embedding dimension %D", dE);
2039: #endif
2040:   /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
2041:   if (dim == dE) {
2042:     for (v = 0; v < Nv; ++v) {
2043:       for (c = 0; c < Nc; ++c) {
2044:         switch (dim)
2045:         {
2046:           case 1: vals[(v*Nc+c)*dim*dim] *= PetscSqr(fegeom->invJ[0]);break;
2047:           case 2: DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break;
2048:           case 3: DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v*Nc+c)*dim*dim], &vals[(v*Nc+c)*dim*dim]);break;
2049:           default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %D for transformation", dim);
2050:         }
2051:       }
2052:     }
2053:   } else {
2054:     for (v = 0; v < Nv; ++v) {
2055:       for (c = 0; c < Nc; ++c) {
2056:         DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v*Nc+c)*dE*dE], &vals[(v*Nc+c)*dE*dE]);
2057:       }
2058:     }
2059:   }
2060:   /* Assume its a vector, otherwise assume its a bunch of scalars */
2061:   if (Nc == 1 || Nc != dim) return(0);
2062:   switch (trans) {
2063:     case IDENTITY_TRANSFORM: break;
2064:     case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
2065:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2066:     case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2067:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2068:   }
2069:   return(0);
2070: }

2072: /*@C
2073:   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.

2075:   Input Parameters:
2076: + dsp        - The PetscDualSpace
2077: . fegeom     - The geometry for this cell
2078: . Nq         - The number of function samples
2079: . Nc         - The number of function components
2080: - pointEval  - The function values

2082:   Output Parameter:
2083: . pointEval  - The transformed function values

2085:   Level: advanced

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

2089:   Note: This only handles tranformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

2091: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2092: @*/
2093: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2094: {
2095:   PetscDualSpaceTransformType trans;
2096:   PetscInt                    k;
2097:   PetscErrorCode              ierr;

2103:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2104:      This determines their transformation properties. */
2105:   PetscDualSpaceGetDeRahm(dsp, &k);
2106:   switch (k)
2107:   {
2108:     case 0: /* H^1 point evaluations */
2109:     trans = IDENTITY_TRANSFORM;break;
2110:     case 1: /* Hcurl preserves tangential edge traces  */
2111:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2112:     case 2:
2113:     case 3: /* Hdiv preserve normal traces */
2114:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2115:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2116:   }
2117:   PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2118:   return(0);
2119: }

2121: /*@C
2122:   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.

2124:   Input Parameters:
2125: + dsp        - The PetscDualSpace
2126: . fegeom     - The geometry for this cell
2127: . Nq         - The number of function samples
2128: . Nc         - The number of function components
2129: - pointEval  - The function values

2131:   Output Parameter:
2132: . pointEval  - The transformed function values

2134:   Level: advanced

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

2138:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

2140: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2141: @*/
2142: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2143: {
2144:   PetscDualSpaceTransformType trans;
2145:   PetscInt                    k;
2146:   PetscErrorCode              ierr;

2152:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2153:      This determines their transformation properties. */
2154:   PetscDualSpaceGetDeRahm(dsp, &k);
2155:   switch (k)
2156:   {
2157:     case 0: /* H^1 point evaluations */
2158:     trans = IDENTITY_TRANSFORM;break;
2159:     case 1: /* Hcurl preserves tangential edge traces  */
2160:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2161:     case 2:
2162:     case 3: /* Hdiv preserve normal traces */
2163:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2164:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2165:   }
2166:   PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2167:   return(0);
2168: }

2170: /*@C
2171:   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.

2173:   Input Parameters:
2174: + dsp        - The PetscDualSpace
2175: . fegeom     - The geometry for this cell
2176: . Nq         - The number of function gradient samples
2177: . Nc         - The number of function components
2178: - pointEval  - The function gradient values

2180:   Output Parameter:
2181: . pointEval  - The transformed function gradient values

2183:   Level: advanced

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

2187:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

2189: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2190: @*/
2191: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2192: {
2193:   PetscDualSpaceTransformType trans;
2194:   PetscInt                    k;
2195:   PetscErrorCode              ierr;

2201:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2202:      This determines their transformation properties. */
2203:   PetscDualSpaceGetDeRahm(dsp, &k);
2204:   switch (k)
2205:   {
2206:     case 0: /* H^1 point evaluations */
2207:     trans = IDENTITY_TRANSFORM;break;
2208:     case 1: /* Hcurl preserves tangential edge traces  */
2209:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2210:     case 2:
2211:     case 3: /* Hdiv preserve normal traces */
2212:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2213:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2214:   }
2215:   PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2216:   return(0);
2217: }

2219: /*@C
2220:   PetscDualSpacePushforwardHessian - Transform the given function Hessian 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.

2222:   Input Parameters:
2223: + dsp        - The PetscDualSpace
2224: . fegeom     - The geometry for this cell
2225: . Nq         - The number of function Hessian samples
2226: . Nc         - The number of function components
2227: - pointEval  - The function gradient values

2229:   Output Parameter:
2230: . pointEval  - The transformed function Hessian values

2232:   Level: advanced

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

2236:   Note: This only handles transformations when the embedding dimension of the geometry in fegeom is the same as the reference dimension.

2238: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2239: @*/
2240: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2241: {
2242:   PetscDualSpaceTransformType trans;
2243:   PetscInt                    k;
2244:   PetscErrorCode              ierr;

2250:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2251:      This determines their transformation properties. */
2252:   PetscDualSpaceGetDeRahm(dsp, &k);
2253:   switch (k)
2254:   {
2255:     case 0: /* H^1 point evaluations */
2256:     trans = IDENTITY_TRANSFORM;break;
2257:     case 1: /* Hcurl preserves tangential edge traces  */
2258:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2259:     case 2:
2260:     case 3: /* Hdiv preserve normal traces */
2261:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2262:     default: SETERRQ1(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2263:   }
2264:   PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2265:   return(0);
2266: }