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: /*
 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}, {0,2}.

 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: `PetscDualSpaceType`, `PetscDualSpaceTensorPointLexicographic_Internal()`
 27: */
 28: PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
 29: {
 30:   while (len--) {
 31:     max -= tup[len];
 32:     if (!max) {
 33:       tup[len] = 0;
 34:       break;
 35:     }
 36:   }
 37:   tup[++len]++;
 38:   return 0;
 39: }

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

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

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

 54:   Level: developer

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

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

 73: /*@C
 74:   PetscDualSpaceRegister - Adds a new `PetscDualSpaceType`

 76:   Not Collective

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

 82:   Sample usage:
 83: .vb
 84:     PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
 85: .ve

 87:   Then, your PetscDualSpace type can be chosen with the procedural interface via
 88: .vb
 89:     PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
 90:     PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
 91: .ve
 92:    or at runtime via the option
 93: .vb
 94:     -petscdualspace_type my_dual_space
 95: .ve

 97:   Level: advanced

 99:   Note:
100:   `PetscDualSpaceRegister()` may be called multiple times to add several user-defined `PetscDualSpace`

102: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceRegisterAll()`, `PetscDualSpaceRegisterDestroy()`
103: @*/
104: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
105: {
106:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
107:   return 0;
108: }

110: /*@C
111:   PetscDualSpaceSetType - Builds a particular `PetscDualSpace` based on its `PetscDualSpaceType`

113:   Collective on sp

115:   Input Parameters:
116: + sp   - The `PetscDualSpace` object
117: - name - The kind of space

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

122:   Level: intermediate

124: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceGetType()`, `PetscDualSpaceCreate()`
125: @*/
126: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
127: {
128:   PetscErrorCode (*r)(PetscDualSpace);
129:   PetscBool match;

132:   PetscObjectTypeCompare((PetscObject)sp, name, &match);
133:   if (match) return 0;

135:   if (!PetscDualSpaceRegisterAllCalled) PetscDualSpaceRegisterAll();
136:   PetscFunctionListFind(PetscDualSpaceList, name, &r);

139:   PetscTryTypeMethod(sp, destroy);
140:   sp->ops->destroy = NULL;

142:   (*r)(sp);
143:   PetscObjectChangeTypeName((PetscObject)sp, name);
144:   return 0;
145: }

147: /*@C
148:   PetscDualSpaceGetType - Gets the `PetscDualSpaceType` name (as a string) from the object.

150:   Not Collective

152:   Input Parameter:
153: . sp  - The `PetscDualSpace`

155:   Output Parameter:
156: . name - The `PetscDualSpaceType` name

158:   Level: intermediate

160: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceSetType()`, `PetscDualSpaceCreate()`
161: @*/
162: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
163: {
166:   if (!PetscDualSpaceRegisterAllCalled) PetscDualSpaceRegisterAll();
167:   *name = ((PetscObject)sp)->type_name;
168:   return 0;
169: }

171: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
172: {
173:   PetscViewerFormat format;
174:   PetscInt          pdim, f;

176:   PetscDualSpaceGetDimension(sp, &pdim);
177:   PetscObjectPrintClassNamePrefixType((PetscObject)sp, v);
178:   PetscViewerASCIIPushTab(v);
179:   if (sp->k) {
180:     PetscViewerASCIIPrintf(v, "Dual space for %" PetscInt_FMT "-forms %swith %" PetscInt_FMT " components, size %" PetscInt_FMT "\n", PetscAbsInt(sp->k), sp->k < 0 ? "(stored in dual form) " : "", sp->Nc, pdim);
181:   } else {
182:     PetscViewerASCIIPrintf(v, "Dual space with %" PetscInt_FMT " components, size %" PetscInt_FMT "\n", sp->Nc, pdim);
183:   }
184:   PetscTryTypeMethod(sp, view, v);
185:   PetscViewerGetFormat(v, &format);
186:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
187:     PetscViewerASCIIPushTab(v);
188:     for (f = 0; f < pdim; ++f) {
189:       PetscViewerASCIIPrintf(v, "Dual basis vector %" PetscInt_FMT "\n", f);
190:       PetscViewerASCIIPushTab(v);
191:       PetscQuadratureView(sp->functional[f], v);
192:       PetscViewerASCIIPopTab(v);
193:     }
194:     PetscViewerASCIIPopTab(v);
195:   }
196:   PetscViewerASCIIPopTab(v);
197:   return 0;
198: }

200: /*@C
201:    PetscDualSpaceViewFromOptions - View a `PetscDualSpace` based on values in the options database

203:    Collective on A

205:    Input Parameters:
206: +  A - the `PetscDualSpace` object
207: .  obj - Optional object, provides the options prefix
208: -  name - command line option name

210:    Level: intermediate

212: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscObjectViewFromOptions()`, `PetscDualSpaceCreate()`
213: @*/
214: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace A, PetscObject obj, const char name[])
215: {
217:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
218:   return 0;
219: }

221: /*@
222:   PetscDualSpaceView - Views a `PetscDualSpace`

224:   Collective on sp

226:   Input Parameters:
227: + sp - the `PetscDualSpace` object to view
228: - v  - the viewer

230:   Level: beginner

232: .seealso: `PetscViewer`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
233: @*/
234: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
235: {
236:   PetscBool iascii;

240:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &v);
241:   PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii);
242:   if (iascii) PetscDualSpaceView_ASCII(sp, v);
243:   return 0;
244: }

246: /*@
247:   PetscDualSpaceSetFromOptions - sets parameters in a `PetscDualSpace` from the options database

249:   Collective on sp

251:   Input Parameter:
252: . sp - the `PetscDualSpace` object to set options for

254:   Options Database Keys:
255: + -petscdualspace_order <order>      - the approximation order of the space
256: . -petscdualspace_form_degree <deg>  - the form degree, say 0 for point evaluations, or 2 for area integrals
257: . -petscdualspace_components <c>     - the number of components, say d for a vector field
258: - -petscdualspace_refcell <celltype> - Reference cell type name

260:   Level: intermediate

262: .seealso: `PetscDualSpaceView()`, `PetscDualSpace`, `PetscObjectSetFromOptions()`
263: @*/
264: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
265: {
266:   DMPolytopeType refCell = DM_POLYTOPE_TRIANGLE;
267:   const char    *defaultType;
268:   char           name[256];
269:   PetscBool      flg;

272:   if (!((PetscObject)sp)->type_name) {
273:     defaultType = PETSCDUALSPACELAGRANGE;
274:   } else {
275:     defaultType = ((PetscObject)sp)->type_name;
276:   }
277:   if (!PetscSpaceRegisterAllCalled) PetscSpaceRegisterAll();

279:   PetscObjectOptionsBegin((PetscObject)sp);
280:   PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
281:   if (flg) {
282:     PetscDualSpaceSetType(sp, name);
283:   } else if (!((PetscObject)sp)->type_name) {
284:     PetscDualSpaceSetType(sp, defaultType);
285:   }
286:   PetscOptionsBoundedInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL, 0);
287:   PetscOptionsInt("-petscdualspace_form_degree", "The form degree of the dofs", "PetscDualSpaceSetFormDegree", sp->k, &sp->k, NULL);
288:   PetscOptionsBoundedInt("-petscdualspace_components", "The number of components", "PetscDualSpaceSetNumComponents", sp->Nc, &sp->Nc, NULL, 1);
289:   PetscTryTypeMethod(sp, setfromoptions, PetscOptionsObject);
290:   PetscOptionsEnum("-petscdualspace_refcell", "Reference cell shape", "PetscDualSpaceSetReferenceCell", DMPolytopeTypes, (PetscEnum)refCell, (PetscEnum *)&refCell, &flg);
291:   if (flg) {
292:     DM K;

294:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K);
295:     PetscDualSpaceSetDM(sp, K);
296:     DMDestroy(&K);
297:   }

299:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
300:   PetscObjectProcessOptionsHandlers((PetscObject)sp, PetscOptionsObject);
301:   PetscOptionsEnd();
302:   sp->setfromoptionscalled = PETSC_TRUE;
303:   return 0;
304: }

306: /*@
307:   PetscDualSpaceSetUp - Construct a basis for a `PetscDualSpace`

309:   Collective on sp

311:   Input Parameter:
312: . sp - the `PetscDualSpace` object to setup

314:   Level: intermediate

316: .seealso: `PetscDualSpaceView()`, `PetscDualSpaceDestroy()`, `PetscDualSpace`
317: @*/
318: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
319: {
321:   if (sp->setupcalled) return 0;
322:   PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
323:   sp->setupcalled = PETSC_TRUE;
324:   PetscTryTypeMethod(sp, setup);
325:   PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
326:   if (sp->setfromoptionscalled) PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
327:   return 0;
328: }

330: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
331: {
332:   PetscInt pStart = -1, pEnd = -1, depth = -1;

334:   if (!dm) return 0;
335:   DMPlexGetChart(dm, &pStart, &pEnd);
336:   DMPlexGetDepth(dm, &depth);

338:   if (sp->pointSpaces) {
339:     PetscInt i;

341:     for (i = 0; i < pEnd - pStart; i++) PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
342:   }
343:   PetscFree(sp->pointSpaces);

345:   if (sp->heightSpaces) {
346:     PetscInt i;

348:     for (i = 0; i <= depth; i++) PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
349:   }
350:   PetscFree(sp->heightSpaces);

352:   PetscSectionDestroy(&(sp->pointSection));
353:   PetscQuadratureDestroy(&(sp->intNodes));
354:   VecDestroy(&(sp->intDofValues));
355:   VecDestroy(&(sp->intNodeValues));
356:   MatDestroy(&(sp->intMat));
357:   PetscQuadratureDestroy(&(sp->allNodes));
358:   VecDestroy(&(sp->allDofValues));
359:   VecDestroy(&(sp->allNodeValues));
360:   MatDestroy(&(sp->allMat));
361:   PetscFree(sp->numDof);
362:   return 0;
363: }

365: /*@
366:   PetscDualSpaceDestroy - Destroys a `PetscDualSpace` object

368:   Collective on sp

370:   Input Parameter:
371: . sp - the `PetscDualSpace` object to destroy

373:   Level: beginner

375: .seealso: `PetscDualSpace`, `PetscDualSpaceView()`, `PetscDualSpace()`, `PetscDualSpaceCreate()`
376: @*/
377: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
378: {
379:   PetscInt dim, f;
380:   DM       dm;

382:   if (!*sp) return 0;

385:   if (--((PetscObject)(*sp))->refct > 0) {
386:     *sp = NULL;
387:     return 0;
388:   }
389:   ((PetscObject)(*sp))->refct = 0;

391:   PetscDualSpaceGetDimension(*sp, &dim);
392:   dm = (*sp)->dm;

394:   PetscTryTypeMethod((*sp), destroy);
395:   PetscDualSpaceClearDMData_Internal(*sp, dm);

397:   for (f = 0; f < dim; ++f) PetscQuadratureDestroy(&(*sp)->functional[f]);
398:   PetscFree((*sp)->functional);
399:   DMDestroy(&(*sp)->dm);
400:   PetscHeaderDestroy(sp);
401:   return 0;
402: }

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

407:   Collective

409:   Input Parameter:
410: . comm - The communicator for the `PetscDualSpace` object

412:   Output Parameter:
413: . sp - The `PetscDualSpace` object

415:   Level: beginner

417: .seealso: `PetscDualSpace`, `PetscDualSpaceSetType()`, `PETSCDUALSPACELAGRANGE`
418: @*/
419: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
420: {
421:   PetscDualSpace s;

424:   PetscCitationsRegister(FECitation, &FEcite);
425:   *sp = NULL;
426:   PetscFEInitializePackage();

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

430:   s->order       = 0;
431:   s->Nc          = 1;
432:   s->k           = 0;
433:   s->spdim       = -1;
434:   s->spintdim    = -1;
435:   s->uniform     = PETSC_TRUE;
436:   s->setupcalled = PETSC_FALSE;

438:   *sp = s;
439:   return 0;
440: }

442: /*@
443:   PetscDualSpaceDuplicate - Creates a duplicate `PetscDualSpace` object that is not setup.

445:   Collective on sp

447:   Input Parameter:
448: . sp - The original `PetscDualSpace`

450:   Output Parameter:
451: . spNew - The duplicate `PetscDualSpace`

453:   Level: beginner

455: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
456: @*/
457: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
458: {
459:   DM                 dm;
460:   PetscDualSpaceType type;
461:   const char        *name;

465:   PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
466:   PetscObjectGetName((PetscObject)sp, &name);
467:   PetscObjectSetName((PetscObject)*spNew, name);
468:   PetscDualSpaceGetType(sp, &type);
469:   PetscDualSpaceSetType(*spNew, type);
470:   PetscDualSpaceGetDM(sp, &dm);
471:   PetscDualSpaceSetDM(*spNew, dm);

473:   (*spNew)->order   = sp->order;
474:   (*spNew)->k       = sp->k;
475:   (*spNew)->Nc      = sp->Nc;
476:   (*spNew)->uniform = sp->uniform;
477:   PetscTryTypeMethod(sp, duplicate, *spNew);
478:   return 0;
479: }

481: /*@
482:   PetscDualSpaceGetDM - Get the `DM` representing the reference cell of a `PetscDualSpace`

484:   Not collective

486:   Input Parameter:
487: . sp - The `PetscDualSpace`

489:   Output Parameter:
490: . dm - The reference cell, that is a `DM` that consists of a single cell

492:   Level: intermediate

494: .seealso: `PetscDualSpace`, `PetscDualSpaceSetDM()`, `PetscDualSpaceCreate()`
495: @*/
496: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
497: {
500:   *dm = sp->dm;
501:   return 0;
502: }

504: /*@
505:   PetscDualSpaceSetDM - Get the `DM` representing the reference cell

507:   Not collective

509:   Input Parameters:
510: + sp - The `PetscDual`Space
511: - dm - The reference cell

513:   Level: intermediate

515: .seealso: `PetscDualSpace`, `DM`, `PetscDualSpaceGetDM()`, `PetscDualSpaceCreate()`
516: @*/
517: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
518: {
522:   PetscObjectReference((PetscObject)dm);
523:   if (sp->dm && sp->dm != dm) PetscDualSpaceClearDMData_Internal(sp, sp->dm);
524:   DMDestroy(&sp->dm);
525:   sp->dm = dm;
526:   return 0;
527: }

529: /*@
530:   PetscDualSpaceGetOrder - Get the order of the dual space

532:   Not collective

534:   Input Parameter:
535: . sp - The `PetscDualSpace`

537:   Output Parameter:
538: . order - The order

540:   Level: intermediate

542: .seealso: `PetscDualSpace`, `PetscDualSpaceSetOrder()`, `PetscDualSpaceCreate()`
543: @*/
544: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
545: {
548:   *order = sp->order;
549:   return 0;
550: }

552: /*@
553:   PetscDualSpaceSetOrder - Set the order of the dual space

555:   Not collective

557:   Input Parameters:
558: + sp - The `PetscDualSpace`
559: - order - The order

561:   Level: intermediate

563: .seealso: `PetscDualSpace`, `PetscDualSpaceGetOrder()`, `PetscDualSpaceCreate()`
564: @*/
565: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
566: {
569:   sp->order = order;
570:   return 0;
571: }

573: /*@
574:   PetscDualSpaceGetNumComponents - Return the number of components for this space

576:   Input Parameter:
577: . sp - The `PetscDualSpace`

579:   Output Parameter:
580: . Nc - The number of components

582:   Level: intermediate

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

587: .seealso: `PetscDualSpaceSetNumComponents()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
588: @*/
589: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
590: {
593:   *Nc = sp->Nc;
594:   return 0;
595: }

597: /*@
598:   PetscDualSpaceSetNumComponents - Set the number of components for this space

600:   Input Parameters:
601: + sp - The `PetscDualSpace`
602: - order - The number of components

604:   Level: intermediate

606: .seealso: `PetscDualSpaceGetNumComponents()`, `PetscDualSpaceCreate()`, `PetscDualSpace`
607: @*/
608: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
609: {
612:   sp->Nc = Nc;
613:   return 0;
614: }

616: /*@
617:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

619:   Not collective

621:   Input Parameters:
622: + sp - The `PetscDualSpace`
623: - i  - The basis number

625:   Output Parameter:
626: . functional - The basis functional

628:   Level: intermediate

630: .seealso: `PetscDualSpace`, `PetscQuadrature`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceCreate()`
631: @*/
632: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
633: {
634:   PetscInt dim;

638:   PetscDualSpaceGetDimension(sp, &dim);
640:   *functional = sp->functional[i];
641:   return 0;
642: }

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

647:   Not collective

649:   Input Parameter:
650: . sp - The `PetscDualSpace`

652:   Output Parameter:
653: . dim - The dimension

655:   Level: intermediate

657: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
658: @*/
659: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
660: {
663:   if (sp->spdim < 0) {
664:     PetscSection section;

666:     PetscDualSpaceGetSection(sp, &section);
667:     if (section) {
668:       PetscSectionGetStorageSize(section, &(sp->spdim));
669:     } else sp->spdim = 0;
670:   }
671:   *dim = sp->spdim;
672:   return 0;
673: }

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

678:   Not collective

680:   Input Parameter:
681: . sp - The `PetscDualSpace`

683:   Output Parameter:
684: . dim - The dimension

686:   Level: intermediate

688: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
689: @*/
690: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
691: {
694:   if (sp->spintdim < 0) {
695:     PetscSection section;

697:     PetscDualSpaceGetSection(sp, &section);
698:     if (section) {
699:       PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
700:     } else sp->spintdim = 0;
701:   }
702:   *intdim = sp->spintdim;
703:   return 0;
704: }

706: /*@
707:    PetscDualSpaceGetUniform - Whether this dual space is uniform

709:    Not collective

711:    Input Parameters:
712: .  sp - A dual space

714:    Output Parameters:
715: .  uniform - `PETSC_TRUE` if (a) the dual space is the same for each point in a stratum of the reference `DMPLEX`, and
716:              (b) every symmetry of each point in the reference `DMPLEX` is also a symmetry of the point's dual space.

718:    Level: advanced

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

724: .seealso: `PetscDualSpace`, `PetscDualSpaceGetPointSubspace()`, `PetscDualSpaceGetSymmetries()`
725: @*/
726: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
727: {
730:   *uniform = sp->uniform;
731:   return 0;
732: }

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

737:   Not collective

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

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

745:   Level: intermediate

747: .seealso: `PetscDualSpace`, `PetscDualSpaceGetFunctional()`, `PetscDualSpaceCreate()`
748: @*/
749: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
750: {
754:   if (!sp->numDof) {
755:     DM           dm;
756:     PetscInt     depth, d;
757:     PetscSection section;

759:     PetscDualSpaceGetDM(sp, &dm);
760:     DMPlexGetDepth(dm, &depth);
761:     PetscCalloc1(depth + 1, &(sp->numDof));
762:     PetscDualSpaceGetSection(sp, &section);
763:     for (d = 0; d <= depth; d++) {
764:       PetscInt dStart, dEnd;

766:       DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
767:       if (dEnd <= dStart) continue;
768:       PetscSectionGetDof(section, dStart, &(sp->numDof[d]));
769:     }
770:   }
771:   *numDof = sp->numDof;
773:   return 0;
774: }

776: /* create the section of the right size and set a permutation for topological ordering */
777: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
778: {
779:   DM           dm;
780:   PetscInt     pStart, pEnd, cStart, cEnd, c, depth, count, i;
781:   PetscInt    *seen, *perm;
782:   PetscSection section;

784:   dm = sp->dm;
785:   PetscSectionCreate(PETSC_COMM_SELF, &section);
786:   DMPlexGetChart(dm, &pStart, &pEnd);
787:   PetscSectionSetChart(section, pStart, pEnd);
788:   PetscCalloc1(pEnd - pStart, &seen);
789:   PetscMalloc1(pEnd - pStart, &perm);
790:   DMPlexGetDepth(dm, &depth);
791:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
792:   for (c = cStart, count = 0; c < cEnd; c++) {
793:     PetscInt  closureSize = -1, e;
794:     PetscInt *closure     = NULL;

796:     perm[count++]    = c;
797:     seen[c - pStart] = 1;
798:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
799:     for (e = 0; e < closureSize; e++) {
800:       PetscInt point = closure[2 * e];

802:       if (seen[point - pStart]) continue;
803:       perm[count++]        = point;
804:       seen[point - pStart] = 1;
805:     }
806:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
807:   }
809:   for (i = 0; i < pEnd - pStart; i++)
810:     if (perm[i] != i) break;
811:   if (i < pEnd - pStart) {
812:     IS permIS;

814:     ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
815:     ISSetPermutation(permIS);
816:     PetscSectionSetPermutation(section, permIS);
817:     ISDestroy(&permIS);
818:   } else {
819:     PetscFree(perm);
820:   }
821:   PetscFree(seen);
822:   *topSection = section;
823:   return 0;
824: }

826: /* mark boundary points and set up */
827: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
828: {
829:   DM       dm;
830:   DMLabel  boundary;
831:   PetscInt pStart, pEnd, p;

833:   dm = sp->dm;
834:   DMLabelCreate(PETSC_COMM_SELF, "boundary", &boundary);
835:   PetscDualSpaceGetDM(sp, &dm);
836:   DMPlexMarkBoundaryFaces(dm, 1, boundary);
837:   DMPlexLabelComplete(dm, boundary);
838:   DMPlexGetChart(dm, &pStart, &pEnd);
839:   for (p = pStart; p < pEnd; p++) {
840:     PetscInt bval;

842:     DMLabelGetValue(boundary, p, &bval);
843:     if (bval == 1) {
844:       PetscInt dof;

846:       PetscSectionGetDof(section, p, &dof);
847:       PetscSectionSetConstraintDof(section, p, dof);
848:     }
849:   }
850:   DMLabelDestroy(&boundary);
851:   PetscSectionSetUp(section);
852:   return 0;
853: }

855: /*@
856:   PetscDualSpaceGetSection - Create a `PetscSection` over the reference cell with the layout from this space

858:   Collective on sp

860:   Input Parameters:
861: . sp      - The `PetscDualSpace`

863:   Output Parameter:
864: . section - The section

866:   Level: advanced

868: .seealso: `PetscDualSpace`, `PetscSection`, `PetscDualSpaceCreate()`, `DMPLEX`
869: @*/
870: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
871: {
872:   PetscInt pStart, pEnd, p;

874:   if (!sp->dm) {
875:     *section = NULL;
876:     return 0;
877:   }
878:   if (!sp->pointSection) {
879:     /* mark the boundary */
880:     PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
881:     DMPlexGetChart(sp->dm, &pStart, &pEnd);
882:     for (p = pStart; p < pEnd; p++) {
883:       PetscDualSpace psp;

885:       PetscDualSpaceGetPointSubspace(sp, p, &psp);
886:       if (psp) {
887:         PetscInt dof;

889:         PetscDualSpaceGetInteriorDimension(psp, &dof);
890:         PetscSectionSetDof(sp->pointSection, p, dof);
891:       }
892:     }
893:     PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection);
894:   }
895:   *section = sp->pointSection;
896:   return 0;
897: }

899: /* this assumes that all of the point dual spaces store their interior dofs first, which is true when the point DMs
900:  * have one cell */
901: PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace sp, PetscInt sStart, PetscInt sEnd)
902: {
903:   PetscReal   *sv0, *v0, *J;
904:   PetscSection section;
905:   PetscInt     dim, s, k;
906:   DM           dm;

908:   PetscDualSpaceGetDM(sp, &dm);
909:   DMGetDimension(dm, &dim);
910:   PetscDualSpaceGetSection(sp, &section);
911:   PetscMalloc3(dim, &v0, dim, &sv0, dim * dim, &J);
912:   PetscDualSpaceGetFormDegree(sp, &k);
913:   for (s = sStart; s < sEnd; s++) {
914:     PetscReal      detJ, hdetJ;
915:     PetscDualSpace ssp;
916:     PetscInt       dof, off, f, sdim;
917:     PetscInt       i, j;
918:     DM             sdm;

920:     PetscDualSpaceGetPointSubspace(sp, s, &ssp);
921:     if (!ssp) continue;
922:     PetscSectionGetDof(section, s, &dof);
923:     PetscSectionGetOffset(section, s, &off);
924:     /* get the first vertex of the reference cell */
925:     PetscDualSpaceGetDM(ssp, &sdm);
926:     DMGetDimension(sdm, &sdim);
927:     DMPlexComputeCellGeometryAffineFEM(sdm, 0, sv0, NULL, NULL, &hdetJ);
928:     DMPlexComputeCellGeometryAffineFEM(dm, s, v0, J, NULL, &detJ);
929:     /* compactify Jacobian */
930:     for (i = 0; i < dim; i++)
931:       for (j = 0; j < sdim; j++) J[i * sdim + j] = J[i * dim + j];
932:     for (f = 0; f < dof; f++) {
933:       PetscQuadrature fn;

935:       PetscDualSpaceGetFunctional(ssp, f, &fn);
936:       PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off + f]));
937:     }
938:   }
939:   PetscFree3(v0, sv0, J);
940:   return 0;
941: }

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

946:   Input Parameters:
947: + sp      - The `PetscDualSpace` object
948: . f       - The basis functional index
949: . time    - The time
950: . 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)
951: . numComp - The number of components for the function
952: . func    - The input function
953: - ctx     - A context for the function

955:   Output Parameter:
956: . value   - numComp output values

958:   Calling Sequence of func:
959: .vb
960:   func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], void *ctx)
961: .ve

963:   Level: beginner

965: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
966: @*/
967: 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)
968: {
972:   PetscUseTypeMethod(sp, apply, f, time, cgeom, numComp, func, ctx, value);
973:   return 0;
974: }

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

979:   Input Parameters:
980: + sp        - The `PetscDualSpace` object
981: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`

983:   Output Parameter:
984: . spValue   - The values of all dual space functionals

986:   Level: advanced

988: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
989: @*/
990: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
991: {
993:   PetscUseTypeMethod(sp, applyall, pointEval, spValue);
994:   return 0;
995: }

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

1000:   Input Parameters:
1001: + sp        - The `PetscDualSpace` object
1002: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`

1004:   Output Parameter:
1005: . spValue   - The values of interior dual space functionals

1007:   Level: advanced

1009: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1010: @*/
1011: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1012: {
1014:   PetscUseTypeMethod(sp, applyint, pointEval, spValue);
1015:   return 0;
1016: }

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

1021:   Input Parameters:
1022: + sp    - The `PetscDualSpace` object
1023: . f     - The basis functional index
1024: . time  - The time
1025: . cgeom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1026: . Nc    - The number of components for the function
1027: . func  - The input function
1028: - ctx   - A context for the function

1030:   Output Parameter:
1031: . value   - The output value

1033:   Calling Sequence of func:
1034: .vb
1035:    func(PetscInt dim, PetscReal time, const PetscReal x[],PetscInt numComponents, PetscScalar values[], void *ctx)
1036: .ve

1038:   Level: advanced

1040:   Note:
1041:   The idea is to evaluate the functional as an integral $ n(f) = \int dx n(x) . f(x) $ where both n and f have Nc components.

1043: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1044: @*/
1045: 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)
1046: {
1047:   DM               dm;
1048:   PetscQuadrature  n;
1049:   const PetscReal *points, *weights;
1050:   PetscReal        x[3];
1051:   PetscScalar     *val;
1052:   PetscInt         dim, dE, qNc, c, Nq, q;
1053:   PetscBool        isAffine;

1057:   PetscDualSpaceGetDM(sp, &dm);
1058:   PetscDualSpaceGetFunctional(sp, f, &n);
1059:   PetscQuadratureGetData(n, &dim, &qNc, &Nq, &points, &weights);
1062:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1063:   *value   = 0.0;
1064:   isAffine = cgeom->isAffine;
1065:   dE       = cgeom->dimEmbed;
1066:   for (q = 0; q < Nq; ++q) {
1067:     if (isAffine) {
1068:       CoordinatesRefToReal(dE, cgeom->dim, cgeom->xi, cgeom->v, cgeom->J, &points[q * dim], x);
1069:       (*func)(dE, time, x, Nc, val, ctx);
1070:     } else {
1071:       (*func)(dE, time, &cgeom->v[dE * q], Nc, val, ctx);
1072:     }
1073:     for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1074:   }
1075:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1076:   return 0;
1077: }

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

1082:   Input Parameters:
1083: + sp        - The `PetscDualSpace` object
1084: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetAllData()`

1086:   Output Parameter:
1087: . spValue   - The values of all dual space functionals

1089:   Level: advanced

1091: .seealso:  `PetscDualSpace`, `PetscDualSpaceCreate()`
1092: @*/
1093: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1094: {
1095:   Vec pointValues, dofValues;
1096:   Mat allMat;

1101:   PetscDualSpaceGetAllData(sp, NULL, &allMat);
1102:   if (!(sp->allNodeValues)) MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1103:   pointValues = sp->allNodeValues;
1104:   if (!(sp->allDofValues)) MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1105:   dofValues = sp->allDofValues;
1106:   VecPlaceArray(pointValues, pointEval);
1107:   VecPlaceArray(dofValues, spValue);
1108:   MatMult(allMat, pointValues, dofValues);
1109:   VecResetArray(dofValues);
1110:   VecResetArray(pointValues);
1111:   return 0;
1112: }

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

1117:   Input Parameters:
1118: + sp        - The `PetscDualSpace` object
1119: - pointEval - Evaluation at the points returned by `PetscDualSpaceGetInteriorData()`

1121:   Output Parameter:
1122: . spValue   - The values of interior dual space functionals

1124:   Level: advanced

1126: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1127: @*/
1128: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1129: {
1130:   Vec pointValues, dofValues;
1131:   Mat intMat;

1136:   PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1137:   if (!(sp->intNodeValues)) MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1138:   pointValues = sp->intNodeValues;
1139:   if (!(sp->intDofValues)) MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1140:   dofValues = sp->intDofValues;
1141:   VecPlaceArray(pointValues, pointEval);
1142:   VecPlaceArray(dofValues, spValue);
1143:   MatMult(intMat, pointValues, dofValues);
1144:   VecResetArray(dofValues);
1145:   VecResetArray(pointValues);
1146:   return 0;
1147: }

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

1152:   Input Parameter:
1153: . sp - The dualspace

1155:   Output Parameters:
1156: + allNodes - A `PetscQuadrature` object containing all evaluation nodes
1157: - allMat - A `Mat` for the node-to-dof transformation

1159:   Level: advanced

1161: .seealso: `PetscQuadrature`, `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`
1162: @*/
1163: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1164: {
1168:   if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1169:     PetscQuadrature qpoints;
1170:     Mat             amat;

1172:     PetscUseTypeMethod(sp, createalldata, &qpoints, &amat);
1173:     PetscQuadratureDestroy(&(sp->allNodes));
1174:     MatDestroy(&(sp->allMat));
1175:     sp->allNodes = qpoints;
1176:     sp->allMat   = amat;
1177:   }
1178:   if (allNodes) *allNodes = sp->allNodes;
1179:   if (allMat) *allMat = sp->allMat;
1180:   return 0;
1181: }

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

1186:   Input Parameter:
1187: . sp - The dualspace

1189:   Output Parameters:
1190: + allNodes - A `PetscQuadrature` object containing all evaluation nodes
1191: - allMat - A `Mat` for the node-to-dof transformation

1193:   Level: advanced

1195: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`, `Mat`, `PetscQuadrature`
1196: @*/
1197: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1198: {
1199:   PetscInt        spdim;
1200:   PetscInt        numPoints, offset;
1201:   PetscReal      *points;
1202:   PetscInt        f, dim;
1203:   PetscInt        Nc, nrows, ncols;
1204:   PetscInt        maxNumPoints;
1205:   PetscQuadrature q;
1206:   Mat             A;

1208:   PetscDualSpaceGetNumComponents(sp, &Nc);
1209:   PetscDualSpaceGetDimension(sp, &spdim);
1210:   if (!spdim) {
1211:     PetscQuadratureCreate(PETSC_COMM_SELF, allNodes);
1212:     PetscQuadratureSetData(*allNodes, 0, 0, 0, NULL, NULL);
1213:   }
1214:   nrows = spdim;
1215:   PetscDualSpaceGetFunctional(sp, 0, &q);
1216:   PetscQuadratureGetData(q, &dim, NULL, &numPoints, NULL, NULL);
1217:   maxNumPoints = numPoints;
1218:   for (f = 1; f < spdim; f++) {
1219:     PetscInt Np;

1221:     PetscDualSpaceGetFunctional(sp, f, &q);
1222:     PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL);
1223:     numPoints += Np;
1224:     maxNumPoints = PetscMax(maxNumPoints, Np);
1225:   }
1226:   ncols = numPoints * Nc;
1227:   PetscMalloc1(dim * numPoints, &points);
1228:   MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1229:   for (f = 0, offset = 0; f < spdim; f++) {
1230:     const PetscReal *p, *w;
1231:     PetscInt         Np, i;
1232:     PetscInt         fnc;

1234:     PetscDualSpaceGetFunctional(sp, f, &q);
1235:     PetscQuadratureGetData(q, NULL, &fnc, &Np, &p, &w);
1237:     for (i = 0; i < Np * dim; i++) points[offset * dim + i] = p[i];
1238:     for (i = 0; i < Np * Nc; i++) MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1239:     offset += Np;
1240:   }
1241:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1242:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1243:   PetscQuadratureCreate(PETSC_COMM_SELF, allNodes);
1244:   PetscQuadratureSetData(*allNodes, dim, 0, numPoints, points, NULL);
1245:   *allMat = A;
1246:   return 0;
1247: }

1249: /*@
1250:   PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1251:   this space, as well as the matrix that computes the degrees of freedom from the quadrature values.  Degrees of
1252:   freedom are interior degrees of freedom if they belong (by `PetscDualSpaceGetSection()`) to interior points in the
1253:   reference `DMPLEX`: complementary boundary degrees of freedom are marked as constrained in the section returned by
1254:   `PetscDualSpaceGetSection()`).

1256:   Input Parameter:
1257: . sp - The dualspace

1259:   Output Parameters:
1260: + intNodes - A `PetscQuadrature` object containing all evaluation points needed to evaluate interior degrees of freedom
1261: - intMat   - A matrix that computes dual space values from point values: size [spdim0 x (npoints * nc)], where spdim0 is
1262:              the size of the constrained layout (`PetscSectionGetConstrainStorageSize()`) of the dual space section,
1263:              npoints is the number of points in intNodes and nc is `PetscDualSpaceGetNumComponents()`.

1265:   Level: advanced

1267: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetDimension()`, `PetscDualSpaceGetNumComponents()`, `PetscQuadratureGetData()`
1268: @*/
1269: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1270: {
1274:   if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1275:     PetscQuadrature qpoints;
1276:     Mat             imat;

1278:     PetscUseTypeMethod(sp, createintdata, &qpoints, &imat);
1279:     PetscQuadratureDestroy(&(sp->intNodes));
1280:     MatDestroy(&(sp->intMat));
1281:     sp->intNodes = qpoints;
1282:     sp->intMat   = imat;
1283:   }
1284:   if (intNodes) *intNodes = sp->intNodes;
1285:   if (intMat) *intMat = sp->intMat;
1286:   return 0;
1287: }

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

1292:   Input Parameter:
1293: . sp - The dualspace

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

1301:   Level: advanced

1303: .seealso: `PetscDualSpace`, `PetscQuadrature`, `Mat`, `PetscDualSpaceCreate()`, `PetscDualSpaceGetInteriorData()`
1304: @*/
1305: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1306: {
1307:   DM              dm;
1308:   PetscInt        spdim0;
1309:   PetscInt        Nc;
1310:   PetscInt        pStart, pEnd, p, f;
1311:   PetscSection    section;
1312:   PetscInt        numPoints, offset, matoffset;
1313:   PetscReal      *points;
1314:   PetscInt        dim;
1315:   PetscInt       *nnz;
1316:   PetscQuadrature q;
1317:   Mat             imat;

1320:   PetscDualSpaceGetSection(sp, &section);
1321:   PetscSectionGetConstrainedStorageSize(section, &spdim0);
1322:   if (!spdim0) {
1323:     *intNodes = NULL;
1324:     *intMat   = NULL;
1325:     return 0;
1326:   }
1327:   PetscDualSpaceGetNumComponents(sp, &Nc);
1328:   PetscSectionGetChart(section, &pStart, &pEnd);
1329:   PetscDualSpaceGetDM(sp, &dm);
1330:   DMGetDimension(dm, &dim);
1331:   PetscMalloc1(spdim0, &nnz);
1332:   for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1333:     PetscInt dof, cdof, off, d;

1335:     PetscSectionGetDof(section, p, &dof);
1336:     PetscSectionGetConstraintDof(section, p, &cdof);
1337:     if (!(dof - cdof)) continue;
1338:     PetscSectionGetOffset(section, p, &off);
1339:     for (d = 0; d < dof; d++, off++, f++) {
1340:       PetscInt Np;

1342:       PetscDualSpaceGetFunctional(sp, off, &q);
1343:       PetscQuadratureGetData(q, NULL, NULL, &Np, NULL, NULL);
1344:       nnz[f] = Np * Nc;
1345:       numPoints += Np;
1346:     }
1347:   }
1348:   MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1349:   PetscFree(nnz);
1350:   PetscMalloc1(dim * numPoints, &points);
1351:   for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1352:     PetscInt dof, cdof, off, d;

1354:     PetscSectionGetDof(section, p, &dof);
1355:     PetscSectionGetConstraintDof(section, p, &cdof);
1356:     if (!(dof - cdof)) continue;
1357:     PetscSectionGetOffset(section, p, &off);
1358:     for (d = 0; d < dof; d++, off++, f++) {
1359:       const PetscReal *p;
1360:       const PetscReal *w;
1361:       PetscInt         Np, i;

1363:       PetscDualSpaceGetFunctional(sp, off, &q);
1364:       PetscQuadratureGetData(q, NULL, NULL, &Np, &p, &w);
1365:       for (i = 0; i < Np * dim; i++) points[offset + i] = p[i];
1366:       for (i = 0; i < Np * Nc; i++) MatSetValue(imat, f, matoffset + i, w[i], INSERT_VALUES);
1367:       offset += Np * dim;
1368:       matoffset += Np * Nc;
1369:     }
1370:   }
1371:   PetscQuadratureCreate(PETSC_COMM_SELF, intNodes);
1372:   PetscQuadratureSetData(*intNodes, dim, 0, numPoints, points, NULL);
1373:   MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1374:   MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1375:   *intMat = imat;
1376:   return 0;
1377: }

1379: /*@
1380:   PetscDualSpaceEqual - Determine if two dual spaces are equivalent

1382:   Input Parameters:
1383: + A    - A `PetscDualSpace` object
1384: - B    - Another `PetscDualSpace` object

1386:   Output Parameter:
1387: . equal - `PETSC_TRUE` if the dual spaces are equivalent

1389:   Level: advanced

1391: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1392: @*/
1393: PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal)
1394: {
1395:   PetscInt        sizeA, sizeB, dimA, dimB;
1396:   const PetscInt *dofA, *dofB;
1397:   PetscQuadrature quadA, quadB;
1398:   Mat             matA, matB;

1403:   *equal = PETSC_FALSE;
1404:   PetscDualSpaceGetDimension(A, &sizeA);
1405:   PetscDualSpaceGetDimension(B, &sizeB);
1406:   if (sizeB != sizeA) return 0;
1407:   DMGetDimension(A->dm, &dimA);
1408:   DMGetDimension(B->dm, &dimB);
1409:   if (dimA != dimB) return 0;

1411:   PetscDualSpaceGetNumDof(A, &dofA);
1412:   PetscDualSpaceGetNumDof(B, &dofB);
1413:   for (PetscInt d = 0; d < dimA; d++) {
1414:     if (dofA[d] != dofB[d]) return 0;
1415:   }

1417:   PetscDualSpaceGetInteriorData(A, &quadA, &matA);
1418:   PetscDualSpaceGetInteriorData(B, &quadB, &matB);
1419:   if (!quadA && !quadB) {
1420:     *equal = PETSC_TRUE;
1421:   } else if (quadA && quadB) {
1422:     PetscQuadratureEqual(quadA, quadB, equal);
1423:     if (*equal == PETSC_FALSE) return 0;
1424:     if (!matA && !matB) return 0;
1425:     if (matA && matB) MatEqual(matA, matB, equal);
1426:     else *equal = PETSC_FALSE;
1427:   }
1428:   return 0;
1429: }

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

1434:   Input Parameters:
1435: + sp    - The `PetscDualSpace` object
1436: . f     - The basis functional index
1437: . time  - The time
1438: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1439: . Nc    - The number of components for the function
1440: . func  - The input function
1441: - ctx   - A context for the function

1443:   Output Parameter:
1444: . value - The output value (scalar)

1446:   Calling Sequence of func:
1447: .vb
1448:   func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt numComponents, PetscScalar values[], void *ctx)
1449: .ve
1450:   Level: advanced

1452:   Note:
1453:   The idea is to evaluate the functional as an integral $ n(f) = \int dx n(x) . f(x)$ where both n and f have Nc components.

1455: .seealso: `PetscDualSpace`, `PetscDualSpaceCreate()`
1456: @*/
1457: 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)
1458: {
1459:   DM               dm;
1460:   PetscQuadrature  n;
1461:   const PetscReal *points, *weights;
1462:   PetscScalar     *val;
1463:   PetscInt         dimEmbed, qNc, c, Nq, q;

1467:   PetscDualSpaceGetDM(sp, &dm);
1468:   DMGetCoordinateDim(dm, &dimEmbed);
1469:   PetscDualSpaceGetFunctional(sp, f, &n);
1470:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1472:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1473:   *value = 0.;
1474:   for (q = 0; q < Nq; ++q) {
1475:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1476:     for (c = 0; c < Nc; ++c) *value += val[c] * weights[q * Nc + c];
1477:   }
1478:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1479:   return 0;
1480: }

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

1486:   Not collective

1488:   Input Parameters:
1489: + sp - the `PetscDualSpace` object
1490: - height - the height of the mesh point for which the subspace is desired

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

1496:   Level: advanced

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

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

1505: .seealso: `PetscDualSpace`, `PetscSpaceGetHeightSubspace()`, `PetscDualSpace`, `PetscDualSpaceGetPointSubspace()`
1506: @*/
1507: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1508: {
1509:   PetscInt depth = -1, cStart, cEnd;
1510:   DM       dm;

1515:   *subsp = NULL;
1516:   dm     = sp->dm;
1517:   DMPlexGetDepth(dm, &depth);
1519:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1520:   if (height == 0 && cEnd == cStart + 1) {
1521:     *subsp = sp;
1522:     return 0;
1523:   }
1524:   if (!sp->heightSpaces) {
1525:     PetscInt h;
1526:     PetscCalloc1(depth + 1, &(sp->heightSpaces));

1528:     for (h = 0; h <= depth; h++) {
1529:       if (h == 0 && cEnd == cStart + 1) continue;
1530:       if (sp->ops->createheightsubspace) (*sp->ops->createheightsubspace)(sp, height, &(sp->heightSpaces[h]));
1531:       else if (sp->pointSpaces) {
1532:         PetscInt hStart, hEnd;

1534:         DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
1535:         if (hEnd > hStart) {
1536:           const char *name;

1538:           PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1539:           if (sp->pointSpaces[hStart]) {
1540:             PetscObjectGetName((PetscObject)sp, &name);
1541:             PetscObjectSetName((PetscObject)sp->pointSpaces[hStart], name);
1542:           }
1543:           sp->heightSpaces[h] = sp->pointSpaces[hStart];
1544:         }
1545:       }
1546:     }
1547:   }
1548:   *subsp = sp->heightSpaces[height];
1549:   return 0;
1550: }

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

1555:   Not collective

1557:   Input Parameters:
1558: + sp - the `PetscDualSpace` object
1559: - point - the point (in the dual space's DM) for which the subspace is desired

1561:   Output Parameters:
1562:   bdsp - the subspace. The functionals in the subspace are with respect to the intrinsic geometry of the
1563:   point, which will be of lesser dimension if height > 0.

1565:   Level: advanced

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

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

1574: .seealso: `PetscDualSpace`, `PetscDualSpaceGetHeightSubspace()`
1575: @*/
1576: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1577: {
1578:   PetscInt pStart = 0, pEnd = 0, cStart, cEnd;
1579:   DM       dm;

1583:   *bdsp = NULL;
1584:   dm    = sp->dm;
1585:   DMPlexGetChart(dm, &pStart, &pEnd);
1587:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1588:   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 */
1589:     *bdsp = sp;
1590:     return 0;
1591:   }
1592:   if (!sp->pointSpaces) {
1593:     PetscInt p;
1594:     PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));

1596:     for (p = 0; p < pEnd - pStart; p++) {
1597:       if (p + pStart == cStart && cEnd == cStart + 1) continue;
1598:       if (sp->ops->createpointsubspace) (*sp->ops->createpointsubspace)(sp, p + pStart, &(sp->pointSpaces[p]));
1599:       else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1600:         PetscInt dim, depth, height;
1601:         DMLabel  label;

1603:         DMPlexGetDepth(dm, &dim);
1604:         DMPlexGetDepthLabel(dm, &label);
1605:         DMLabelGetValue(label, p + pStart, &depth);
1606:         height = dim - depth;
1607:         PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1608:         PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1609:       }
1610:     }
1611:   }
1612:   *bdsp = sp->pointSpaces[point - pStart];
1613:   return 0;
1614: }

1616: /*@C
1617:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1619:   Not collective

1621:   Input Parameter:
1622: . sp - the `PetscDualSpace` object

1624:   Output Parameters:
1625: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1626: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation

1628:   Level: developer

1630:   Note:
1631:   The permutation and flip arrays are organized in the following way
1632: .vb
1633:   perms[p][ornt][dof # on point] = new local dof #
1634:   flips[p][ornt][dof # on point] = reversal or not
1635: .ve

1637: .seealso: `PetscDualSpace`
1638: @*/
1639: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1640: {
1642:   if (perms) {
1644:     *perms = NULL;
1645:   }
1646:   if (flips) {
1648:     *flips = NULL;
1649:   }
1650:   if (sp->ops->getsymmetries) (sp->ops->getsymmetries)(sp, perms, flips);
1651:   return 0;
1652: }

1654: /*@
1655:   PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1656:   dual space's functionals.

1658:   Input Parameter:
1659: . dsp - The `PetscDualSpace`

1661:   Output Parameter:
1662: . k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1663:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1664:         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).
1665:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1666:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1667:         but are stored as 1-forms.

1669:   Level: developer

1671: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1672: @*/
1673: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1674: {
1678:   *k = dsp->k;
1679:   return 0;
1680: }

1682: /*@
1683:   PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1684:   dual space's functionals.

1686:   Input Parameters:
1687: + dsp - The `PetscDualSpace`
1688: - k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1689:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1690:         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).
1691:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1692:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1693:         but are stored as 1-forms.

1695:   Level: developer

1697: .seealso: `PetscDualSpace`, `PetscDTAltV`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1698: @*/
1699: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1700: {
1701:   PetscInt dim;

1706:   dim = dsp->dm->dim;
1708:   dsp->k = k;
1709:   return 0;
1710: }

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

1715:   Input Parameter:
1716: . dsp - The `PetscDualSpace`

1718:   Output Parameter:
1719: . k   - The simplex dimension

1721:   Level: developer

1723:   Note:
1724:   Currently supported values are
1725: .vb
1726:   0: These are H_1 methods that only transform coordinates
1727:   1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1728:   2: These are the same as 1
1729:   3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)
1730: .ve

1732: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceTransformType`
1733: @*/
1734: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1735: {
1736:   PetscInt dim;

1741:   dim = dsp->dm->dim;
1742:   if (!dsp->k) *k = IDENTITY_TRANSFORM;
1743:   else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1744:   else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1745:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1746:   return 0;
1747: }

1749: /*@C
1750:   PetscDualSpaceTransform - Transform the function values

1752:   Input Parameters:
1753: + dsp       - The `PetscDualSpace`
1754: . trans     - The type of transform
1755: . isInverse - Flag to invert the transform
1756: . fegeom    - The cell geometry
1757: . Nv        - The number of function samples
1758: . Nc        - The number of function components
1759: - vals      - The function values

1761:   Output Parameter:
1762: . vals      - The transformed function values

1764:   Level: intermediate

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

1769: .seealso: `PetscDualSpace`, `PetscDualSpaceTransformGradient()`, `PetscDualSpaceTransformHessian()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1770: @*/
1771: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1772: {
1773:   PetscReal Jstar[9] = {0};
1774:   PetscInt  dim, v, c, Nk;

1780:   /* TODO: not handling dimEmbed != dim right now */
1781:   dim = dsp->dm->dim;
1782:   /* No change needed for 0-forms */
1783:   if (!dsp->k) return 0;
1784:   PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1785:   /* TODO: use fegeom->isAffine */
1786:   PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1787:   for (v = 0; v < Nv; ++v) {
1788:     switch (Nk) {
1789:     case 1:
1790:       for (c = 0; c < Nc; c++) vals[v * Nc + c] *= Jstar[0];
1791:       break;
1792:     case 2:
1793:       for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1794:       break;
1795:     case 3:
1796:       for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v * Nc + c], &vals[v * Nc + c]);
1797:       break;
1798:     default:
1799:       SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %" PetscInt_FMT " for transformation", Nk);
1800:     }
1801:   }
1802:   return 0;
1803: }

1805: /*@C
1806:   PetscDualSpaceTransformGradient - Transform the function gradient values

1808:   Input Parameters:
1809: + dsp       - The `PetscDualSpace`
1810: . trans     - The type of transform
1811: . isInverse - Flag to invert the transform
1812: . fegeom    - The cell geometry
1813: . Nv        - The number of function gradient samples
1814: . Nc        - The number of function components
1815: - vals      - The function gradient values

1817:   Output Parameter:
1818: . vals      - The transformed function gradient values

1820:   Level: intermediate

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

1825: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1826: @*/
1827: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1828: {
1829:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1830:   PetscInt       v, c, d;

1836: #ifdef PETSC_USE_DEBUG
1838: #endif
1839:   /* Transform gradient */
1840:   if (dim == dE) {
1841:     for (v = 0; v < Nv; ++v) {
1842:       for (c = 0; c < Nc; ++c) {
1843:         switch (dim) {
1844:         case 1:
1845:           vals[(v * Nc + c) * dim] *= fegeom->invJ[0];
1846:           break;
1847:         case 2:
1848:           DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1849:           break;
1850:         case 3:
1851:           DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, 1, &vals[(v * Nc + c) * dim], &vals[(v * Nc + c) * dim]);
1852:           break;
1853:         default:
1854:           SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1855:         }
1856:       }
1857:     }
1858:   } else {
1859:     for (v = 0; v < Nv; ++v) {
1860:       for (c = 0; c < Nc; ++c) DMPlex_MultTransposeReal_Internal(fegeom->invJ, dim, dE, 1, &vals[(v * Nc + c) * dE], &vals[(v * Nc + c) * dE]);
1861:     }
1862:   }
1863:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1864:   if (Nc == 1 || Nc != dim) return 0;
1865:   switch (trans) {
1866:   case IDENTITY_TRANSFORM:
1867:     break;
1868:   case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
1869:     if (isInverse) {
1870:       for (v = 0; v < Nv; ++v) {
1871:         for (d = 0; d < dim; ++d) {
1872:           switch (dim) {
1873:           case 2:
1874:             DMPlex_MultTranspose2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1875:             break;
1876:           case 3:
1877:             DMPlex_MultTranspose3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1878:             break;
1879:           default:
1880:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1881:           }
1882:         }
1883:       }
1884:     } else {
1885:       for (v = 0; v < Nv; ++v) {
1886:         for (d = 0; d < dim; ++d) {
1887:           switch (dim) {
1888:           case 2:
1889:             DMPlex_MultTranspose2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1890:             break;
1891:           case 3:
1892:             DMPlex_MultTranspose3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1893:             break;
1894:           default:
1895:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1896:           }
1897:         }
1898:       }
1899:     }
1900:     break;
1901:   case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
1902:     if (isInverse) {
1903:       for (v = 0; v < Nv; ++v) {
1904:         for (d = 0; d < dim; ++d) {
1905:           switch (dim) {
1906:           case 2:
1907:             DMPlex_Mult2DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1908:             break;
1909:           case 3:
1910:             DMPlex_Mult3DReal_Internal(fegeom->invJ, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1911:             break;
1912:           default:
1913:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1914:           }
1915:           for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] *= fegeom->detJ[0];
1916:         }
1917:       }
1918:     } else {
1919:       for (v = 0; v < Nv; ++v) {
1920:         for (d = 0; d < dim; ++d) {
1921:           switch (dim) {
1922:           case 2:
1923:             DMPlex_Mult2DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1924:             break;
1925:           case 3:
1926:             DMPlex_Mult3DReal_Internal(fegeom->J, dim, &vals[v * Nc * dim + d], &vals[v * Nc * dim + d]);
1927:             break;
1928:           default:
1929:             SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1930:           }
1931:           for (c = 0; c < Nc; ++c) vals[(v * Nc + c) * dim + d] /= fegeom->detJ[0];
1932:         }
1933:       }
1934:     }
1935:     break;
1936:   }
1937:   return 0;
1938: }

1940: /*@C
1941:   PetscDualSpaceTransformHessian - Transform the function Hessian values

1943:   Input Parameters:
1944: + dsp       - The `PetscDualSpace`
1945: . trans     - The type of transform
1946: . isInverse - Flag to invert the transform
1947: . fegeom    - The cell geometry
1948: . Nv        - The number of function Hessian samples
1949: . Nc        - The number of function components
1950: - vals      - The function gradient values

1952:   Output Parameter:
1953: . vals      - The transformed function Hessian values

1955:   Level: intermediate

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

1960: .seealso: `PetscDualSpace`, `PetscDualSpaceTransform()`, `PetscDualSpacePullback()`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransformType`
1961: @*/
1962: PetscErrorCode PetscDualSpaceTransformHessian(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1963: {
1964:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1965:   PetscInt       v, c;

1971: #ifdef PETSC_USE_DEBUG
1973: #endif
1974:   /* Transform Hessian: J^{-T}_{ik} J^{-T}_{jl} H(f)_{kl} = J^{-T}_{ik} H(f)_{kl} J^{-1}_{lj} */
1975:   if (dim == dE) {
1976:     for (v = 0; v < Nv; ++v) {
1977:       for (c = 0; c < Nc; ++c) {
1978:         switch (dim) {
1979:         case 1:
1980:           vals[(v * Nc + c) * dim * dim] *= PetscSqr(fegeom->invJ[0]);
1981:           break;
1982:         case 2:
1983:           DMPlex_PTAP2DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
1984:           break;
1985:         case 3:
1986:           DMPlex_PTAP3DReal_Internal(fegeom->invJ, &vals[(v * Nc + c) * dim * dim], &vals[(v * Nc + c) * dim * dim]);
1987:           break;
1988:         default:
1989:           SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dim %" PetscInt_FMT " for transformation", dim);
1990:         }
1991:       }
1992:     }
1993:   } else {
1994:     for (v = 0; v < Nv; ++v) {
1995:       for (c = 0; c < Nc; ++c) DMPlex_PTAPReal_Internal(fegeom->invJ, dim, dE, &vals[(v * Nc + c) * dE * dE], &vals[(v * Nc + c) * dE * dE]);
1996:     }
1997:   }
1998:   /* Assume its a vector, otherwise assume its a bunch of scalars */
1999:   if (Nc == 1 || Nc != dim) return 0;
2000:   switch (trans) {
2001:   case IDENTITY_TRANSFORM:
2002:     break;
2003:   case COVARIANT_PIOLA_TRANSFORM: /* Covariant Piola mapping $\sigma^*(F) = J^{-T} F \circ \phi^{-1)$ */
2004:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2005:   case CONTRAVARIANT_PIOLA_TRANSFORM: /* Contravariant Piola mapping $\sigma^*(F) = \frac{1}{|\det J|} J F \circ \phi^{-1}$ */
2006:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Piola mapping for Hessians not yet supported");
2007:   }
2008:   return 0;
2009: }

2011: /*@C
2012:   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.

2014:   Input Parameters:
2015: + dsp        - The `PetscDualSpace`
2016: . fegeom     - The geometry for this cell
2017: . Nq         - The number of function samples
2018: . Nc         - The number of function components
2019: - pointEval  - The function values

2021:   Output Parameter:
2022: . pointEval  - The transformed function values

2024:   Level: advanced

2026:   Notes:
2027:   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.

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

2031: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2032: @*/
2033: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2034: {
2035:   PetscDualSpaceTransformType trans;
2036:   PetscInt                    k;

2042:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2043:      This determines their transformation properties. */
2044:   PetscDualSpaceGetDeRahm(dsp, &k);
2045:   switch (k) {
2046:   case 0: /* H^1 point evaluations */
2047:     trans = IDENTITY_TRANSFORM;
2048:     break;
2049:   case 1: /* Hcurl preserves tangential edge traces  */
2050:     trans = COVARIANT_PIOLA_TRANSFORM;
2051:     break;
2052:   case 2:
2053:   case 3: /* Hdiv preserve normal traces */
2054:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;
2055:     break;
2056:   default:
2057:     SETERRQ(PetscObjectComm((PetscObject)dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %" PetscInt_FMT " for transformation", k);
2058:   }
2059:   PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2060:   return 0;
2061: }

2063: /*@C
2064:   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.

2066:   Input Parameters:
2067: + dsp        - The `PetscDualSpace`
2068: . fegeom     - The geometry for this cell
2069: . Nq         - The number of function samples
2070: . Nc         - The number of function components
2071: - pointEval  - The function values

2073:   Output Parameter:
2074: . pointEval  - The transformed function values

2076:   Level: advanced

2078:   Notes:
2079:   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.

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

2083: .seealso: `PetscDualSpace`, `PetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2084: @*/
2085: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2086: {
2087:   PetscDualSpaceTransformType trans;
2088:   PetscInt                    k;

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

2115: /*@C
2116:   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.

2118:   Input Parameters:
2119: + dsp        - The `PetscDualSpace`
2120: . fegeom     - The geometry for this cell
2121: . Nq         - The number of function gradient samples
2122: . Nc         - The number of function components
2123: - pointEval  - The function gradient values

2125:   Output Parameter:
2126: . pointEval  - The transformed function gradient values

2128:   Level: advanced

2130:   Notes:
2131:   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.

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

2135: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2136: @*/
2137: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2138: {
2139:   PetscDualSpaceTransformType trans;
2140:   PetscInt                    k;

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

2167: /*@C
2168:   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.

2170:   Input Parameters:
2171: + dsp        - The `PetscDualSpace`
2172: . fegeom     - The geometry for this cell
2173: . Nq         - The number of function Hessian samples
2174: . Nc         - The number of function components
2175: - pointEval  - The function gradient values

2177:   Output Parameter:
2178: . pointEval  - The transformed function Hessian values

2180:   Level: advanced

2182:   Notes:
2183:   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.

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

2187: .seealso: `PetscDualSpace`, `PetscDualSpacePushforward()`, `PPetscDualSpacePullback()`, `PetscDualSpaceTransform()`, `PetscDualSpaceGetDeRahm()`
2188: @*/
2189: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2190: {
2191:   PetscDualSpaceTransformType trans;
2192:   PetscInt                    k;

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