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: 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: 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 PetscDualSpace implementation

 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:   Notes:
 83:   PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces

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

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

100:   Level: advanced

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

104: @*/
105: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
106: {
107:   PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
108:   return 0;
109: }

111: /*@C
112:   PetscDualSpaceSetType - Builds a particular PetscDualSpace

114:   Collective on sp

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

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

123:   Level: intermediate

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

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

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

140:   if (sp->ops->destroy) {
141:     (*sp->ops->destroy)(sp);
142:     sp->ops->destroy = NULL;
143:   }
144:   (*r)(sp);
145:   PetscObjectChangeTypeName((PetscObject) sp, name);
146:   return 0;
147: }

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

152:   Not Collective

154:   Input Parameter:
155: . sp  - The PetscDualSpace

157:   Output Parameter:
158: . name - The PetscDualSpace type name

160:   Level: intermediate

162: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
163: @*/
164: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
165: {
168:   if (!PetscDualSpaceRegisterAllCalled) {
169:     PetscDualSpaceRegisterAll();
170:   }
171:   *name = ((PetscObject) sp)->type_name;
172:   return 0;
173: }

175: static PetscErrorCode PetscDualSpaceView_ASCII(PetscDualSpace sp, PetscViewer v)
176: {
177:   PetscViewerFormat format;
178:   PetscInt          pdim, f;

180:   PetscDualSpaceGetDimension(sp, &pdim);
181:   PetscObjectPrintClassNamePrefixType((PetscObject) sp, v);
182:   PetscViewerASCIIPushTab(v);
183:   if (sp->k) {
184:     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);
185:   } else {
186:     PetscViewerASCIIPrintf(v, "Dual space with %D components, size %D\n", sp->Nc, pdim);
187:   }
188:   if (sp->ops->view) (*sp->ops->view)(sp, v);
189:   PetscViewerGetFormat(v, &format);
190:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
191:     PetscViewerASCIIPushTab(v);
192:     for (f = 0; f < pdim; ++f) {
193:       PetscViewerASCIIPrintf(v, "Dual basis vector %D\n", f);
194:       PetscViewerASCIIPushTab(v);
195:       PetscQuadratureView(sp->functional[f], v);
196:       PetscViewerASCIIPopTab(v);
197:     }
198:     PetscViewerASCIIPopTab(v);
199:   }
200:   PetscViewerASCIIPopTab(v);
201:   return 0;
202: }

204: /*@C
205:    PetscDualSpaceViewFromOptions - View from Options

207:    Collective on PetscDualSpace

209:    Input Parameters:
210: +  A - the PetscDualSpace object
211: .  obj - Optional object, proivides prefix
212: -  name - command line option

214:    Level: intermediate
215: .seealso:  PetscDualSpace, PetscDualSpaceView(), PetscObjectViewFromOptions(), PetscDualSpaceCreate()
216: @*/
217: PetscErrorCode  PetscDualSpaceViewFromOptions(PetscDualSpace A,PetscObject obj,const char name[])
218: {
220:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
221:   return 0;
222: }

224: /*@
225:   PetscDualSpaceView - Views a PetscDualSpace

227:   Collective on sp

229:   Input Parameters:
230: + sp - the PetscDualSpace object to view
231: - v  - the viewer

233:   Level: beginner

235: .seealso PetscDualSpaceDestroy(), PetscDualSpace
236: @*/
237: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
238: {
239:   PetscBool      iascii;

243:   if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
244:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
245:   if (iascii) PetscDualSpaceView_ASCII(sp, v);
246:   return 0;
247: }

249: /*@
250:   PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database

252:   Collective on sp

254:   Input Parameter:
255: . sp - the PetscDualSpace object to set options for

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

263:   Level: intermediate

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

276:   if (!((PetscObject) sp)->type_name) {
277:     defaultType = PETSCDUALSPACELAGRANGE;
278:   } else {
279:     defaultType = ((PetscObject) sp)->type_name;
280:   }
281:   if (!PetscSpaceRegisterAllCalled) PetscSpaceRegisterAll();

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

300:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, refCell, &K);
301:     PetscDualSpaceSetDM(sp, K);
302:     DMDestroy(&K);
303:   }

305:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
306:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) sp);
307:   PetscOptionsEnd();
308:   sp->setfromoptionscalled = PETSC_TRUE;
309:   return 0;
310: }

312: /*@
313:   PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace

315:   Collective on sp

317:   Input Parameter:
318: . sp - the PetscDualSpace object to setup

320:   Level: intermediate

322: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy(), PetscDualSpace
323: @*/
324: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
325: {
327:   if (sp->setupcalled) return 0;
328:   PetscLogEventBegin(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
329:   sp->setupcalled = PETSC_TRUE;
330:   if (sp->ops->setup) (*sp->ops->setup)(sp);
331:   PetscLogEventEnd(PETSCDUALSPACE_SetUp, sp, 0, 0, 0);
332:   if (sp->setfromoptionscalled) PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
333:   return 0;
334: }

336: static PetscErrorCode PetscDualSpaceClearDMData_Internal(PetscDualSpace sp, DM dm)
337: {
338:   PetscInt       pStart = -1, pEnd = -1, depth = -1;

340:   if (!dm) return 0;
341:   DMPlexGetChart(dm, &pStart, &pEnd);
342:   DMPlexGetDepth(dm, &depth);

344:   if (sp->pointSpaces) {
345:     PetscInt i;

347:     for (i = 0; i < pEnd - pStart; i++) {
348:       PetscDualSpaceDestroy(&(sp->pointSpaces[i]));
349:     }
350:   }
351:   PetscFree(sp->pointSpaces);

353:   if (sp->heightSpaces) {
354:     PetscInt i;

356:     for (i = 0; i <= depth; i++) {
357:       PetscDualSpaceDestroy(&(sp->heightSpaces[i]));
358:     }
359:   }
360:   PetscFree(sp->heightSpaces);

362:   PetscSectionDestroy(&(sp->pointSection));
363:   PetscQuadratureDestroy(&(sp->intNodes));
364:   VecDestroy(&(sp->intDofValues));
365:   VecDestroy(&(sp->intNodeValues));
366:   MatDestroy(&(sp->intMat));
367:   PetscQuadratureDestroy(&(sp->allNodes));
368:   VecDestroy(&(sp->allDofValues));
369:   VecDestroy(&(sp->allNodeValues));
370:   MatDestroy(&(sp->allMat));
371:   PetscFree(sp->numDof);
372:   return 0;
373: }

375: /*@
376:   PetscDualSpaceDestroy - Destroys a PetscDualSpace object

378:   Collective on sp

380:   Input Parameter:
381: . sp - the PetscDualSpace object to destroy

383:   Level: beginner

385: .seealso PetscDualSpaceView(), PetscDualSpace(), PetscDualSpaceCreate()
386: @*/
387: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
388: {
389:   PetscInt       dim, f;
390:   DM             dm;

392:   if (!*sp) return 0;

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

398:   PetscDualSpaceGetDimension(*sp, &dim);
399:   dm = (*sp)->dm;

401:   if ((*sp)->ops->destroy) (*(*sp)->ops->destroy)(*sp);
402:   PetscDualSpaceClearDMData_Internal(*sp, dm);

404:   for (f = 0; f < dim; ++f) {
405:     PetscQuadratureDestroy(&(*sp)->functional[f]);
406:   }
407:   PetscFree((*sp)->functional);
408:   DMDestroy(&(*sp)->dm);
409:   PetscHeaderDestroy(sp);
410:   return 0;
411: }

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

416:   Collective

418:   Input Parameter:
419: . comm - The communicator for the PetscDualSpace object

421:   Output Parameter:
422: . sp - The PetscDualSpace object

424:   Level: beginner

426: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
427: @*/
428: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
429: {
430:   PetscDualSpace s;

433:   PetscCitationsRegister(FECitation,&FEcite);
434:   *sp  = NULL;
435:   PetscFEInitializePackage();

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

439:   s->order       = 0;
440:   s->Nc          = 1;
441:   s->k           = 0;
442:   s->spdim       = -1;
443:   s->spintdim    = -1;
444:   s->uniform     = PETSC_TRUE;
445:   s->setupcalled = PETSC_FALSE;

447:   *sp = s;
448:   return 0;
449: }

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

454:   Collective on sp

456:   Input Parameter:
457: . sp - The original PetscDualSpace

459:   Output Parameter:
460: . spNew - The duplicate PetscDualSpace

462:   Level: beginner

464: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
465: @*/
466: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
467: {
468:   DM             dm;
469:   PetscDualSpaceType type;
470:   const char     *name;

474:   PetscDualSpaceCreate(PetscObjectComm((PetscObject)sp), spNew);
475:   PetscObjectGetName((PetscObject) sp,     &name);
476:   PetscObjectSetName((PetscObject) *spNew,  name);
477:   PetscDualSpaceGetType(sp, &type);
478:   PetscDualSpaceSetType(*spNew, type);
479:   PetscDualSpaceGetDM(sp, &dm);
480:   PetscDualSpaceSetDM(*spNew, dm);

482:   (*spNew)->order   = sp->order;
483:   (*spNew)->k       = sp->k;
484:   (*spNew)->Nc      = sp->Nc;
485:   (*spNew)->uniform = sp->uniform;
486:   if (sp->ops->duplicate) (*sp->ops->duplicate)(sp, *spNew);
487:   return 0;
488: }

490: /*@
491:   PetscDualSpaceGetDM - Get the DM representing the reference cell

493:   Not collective

495:   Input Parameter:
496: . sp - The PetscDualSpace

498:   Output Parameter:
499: . dm - The reference cell

501:   Level: intermediate

503: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
504: @*/
505: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
506: {
509:   *dm = sp->dm;
510:   return 0;
511: }

513: /*@
514:   PetscDualSpaceSetDM - Get the DM representing the reference cell

516:   Not collective

518:   Input Parameters:
519: + sp - The PetscDualSpace
520: - dm - The reference cell

522:   Level: intermediate

524: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
525: @*/
526: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
527: {
531:   PetscObjectReference((PetscObject) dm);
532:   if (sp->dm && sp->dm != dm) {
533:     PetscDualSpaceClearDMData_Internal(sp, sp->dm);
534:   }
535:   DMDestroy(&sp->dm);
536:   sp->dm = dm;
537:   return 0;
538: }

540: /*@
541:   PetscDualSpaceGetOrder - Get the order of the dual space

543:   Not collective

545:   Input Parameter:
546: . sp - The PetscDualSpace

548:   Output Parameter:
549: . order - The order

551:   Level: intermediate

553: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
554: @*/
555: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
556: {
559:   *order = sp->order;
560:   return 0;
561: }

563: /*@
564:   PetscDualSpaceSetOrder - Set the order of the dual space

566:   Not collective

568:   Input Parameters:
569: + sp - The PetscDualSpace
570: - order - The order

572:   Level: intermediate

574: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
575: @*/
576: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
577: {
580:   sp->order = order;
581:   return 0;
582: }

584: /*@
585:   PetscDualSpaceGetNumComponents - Return the number of components for this space

587:   Input Parameter:
588: . sp - The PetscDualSpace

590:   Output Parameter:
591: . Nc - The number of components

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

595:   Level: intermediate

597: .seealso: PetscDualSpaceSetNumComponents(), PetscDualSpaceGetDimension(), PetscDualSpaceCreate(), PetscDualSpace
598: @*/
599: PetscErrorCode PetscDualSpaceGetNumComponents(PetscDualSpace sp, PetscInt *Nc)
600: {
603:   *Nc = sp->Nc;
604:   return 0;
605: }

607: /*@
608:   PetscDualSpaceSetNumComponents - Set the number of components for this space

610:   Input Parameters:
611: + sp - The PetscDualSpace
612: - order - The number of components

614:   Level: intermediate

616: .seealso: PetscDualSpaceGetNumComponents(), PetscDualSpaceCreate(), PetscDualSpace
617: @*/
618: PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace sp, PetscInt Nc)
619: {
622:   sp->Nc = Nc;
623:   return 0;
624: }

626: /*@
627:   PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space

629:   Not collective

631:   Input Parameters:
632: + sp - The PetscDualSpace
633: - i  - The basis number

635:   Output Parameter:
636: . functional - The basis functional

638:   Level: intermediate

640: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
641: @*/
642: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
643: {
644:   PetscInt       dim;

648:   PetscDualSpaceGetDimension(sp, &dim);
650:   *functional = sp->functional[i];
651:   return 0;
652: }

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

657:   Not collective

659:   Input Parameter:
660: . sp - The PetscDualSpace

662:   Output Parameter:
663: . dim - The dimension

665:   Level: intermediate

667: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
668: @*/
669: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
670: {
673:   if (sp->spdim < 0) {
674:     PetscSection section;

676:     PetscDualSpaceGetSection(sp, &section);
677:     if (section) {
678:       PetscSectionGetStorageSize(section, &(sp->spdim));
679:     } else sp->spdim = 0;
680:   }
681:   *dim = sp->spdim;
682:   return 0;
683: }

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

688:   Not collective

690:   Input Parameter:
691: . sp - The PetscDualSpace

693:   Output Parameter:
694: . dim - The dimension

696:   Level: intermediate

698: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
699: @*/
700: PetscErrorCode PetscDualSpaceGetInteriorDimension(PetscDualSpace sp, PetscInt *intdim)
701: {
704:   if (sp->spintdim < 0) {
705:     PetscSection section;

707:     PetscDualSpaceGetSection(sp, &section);
708:     if (section) {
709:       PetscSectionGetConstrainedStorageSize(section, &(sp->spintdim));
710:     } else sp->spintdim = 0;
711:   }
712:   *intdim = sp->spintdim;
713:   return 0;
714: }

716: /*@
717:    PetscDualSpaceGetUniform - Whether this dual space is uniform

719:    Not collective

721:    Input Parameters:
722: .  sp - A dual space

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

728:    Level: advanced

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

733: .seealso: PetscDualSpaceGetPointSubspace(), PetscDualSpaceGetSymmetries()
734: @*/
735: PetscErrorCode PetscDualSpaceGetUniform(PetscDualSpace sp, PetscBool *uniform)
736: {
739:   *uniform = sp->uniform;
740:   return 0;
741: }

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

746:   Not collective

748:   Input Parameter:
749: . sp - The PetscDualSpace

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

754:   Level: intermediate

756: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
757: @*/
758: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
759: {
763:   if (!sp->numDof) {
764:     DM       dm;
765:     PetscInt depth, d;
766:     PetscSection section;

768:     PetscDualSpaceGetDM(sp, &dm);
769:     DMPlexGetDepth(dm, &depth);
770:     PetscCalloc1(depth+1,&(sp->numDof));
771:     PetscDualSpaceGetSection(sp, &section);
772:     for (d = 0; d <= depth; d++) {
773:       PetscInt dStart, dEnd;

775:       DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
776:       if (dEnd <= dStart) continue;
777:       PetscSectionGetDof(section, dStart, &(sp->numDof[d]));

779:     }
780:   }
781:   *numDof = sp->numDof;
783:   return 0;
784: }

786: /* create the section of the right size and set a permutation for topological ordering */
787: PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace sp, PetscSection *topSection)
788: {
789:   DM             dm;
790:   PetscInt       pStart, pEnd, cStart, cEnd, c, depth, count, i;
791:   PetscInt       *seen, *perm;
792:   PetscSection   section;

794:   dm = sp->dm;
795:   PetscSectionCreate(PETSC_COMM_SELF, &section);
796:   DMPlexGetChart(dm, &pStart, &pEnd);
797:   PetscSectionSetChart(section, pStart, pEnd);
798:   PetscCalloc1(pEnd - pStart, &seen);
799:   PetscMalloc1(pEnd - pStart, &perm);
800:   DMPlexGetDepth(dm, &depth);
801:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
802:   for (c = cStart, count = 0; c < cEnd; c++) {
803:     PetscInt closureSize = -1, e;
804:     PetscInt *closure = NULL;

806:     perm[count++] = c;
807:     seen[c-pStart] = 1;
808:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
809:     for (e = 0; e < closureSize; e++) {
810:       PetscInt point = closure[2*e];

812:       if (seen[point-pStart]) continue;
813:       perm[count++] = point;
814:       seen[point-pStart] = 1;
815:     }
816:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
817:   }
819:   for (i = 0; i < pEnd - pStart; i++) if (perm[i] != i) break;
820:   if (i < pEnd - pStart) {
821:     IS permIS;

823:     ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS);
824:     ISSetPermutation(permIS);
825:     PetscSectionSetPermutation(section, permIS);
826:     ISDestroy(&permIS);
827:   } else {
828:     PetscFree(perm);
829:   }
830:   PetscFree(seen);
831:   *topSection = section;
832:   return 0;
833: }

835: /* mark boundary points and set up */
836: PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace sp, PetscSection section)
837: {
838:   DM             dm;
839:   DMLabel        boundary;
840:   PetscInt       pStart, pEnd, p;

842:   dm = sp->dm;
843:   DMLabelCreate(PETSC_COMM_SELF,"boundary",&boundary);
844:   PetscDualSpaceGetDM(sp,&dm);
845:   DMPlexMarkBoundaryFaces(dm,1,boundary);
846:   DMPlexLabelComplete(dm,boundary);
847:   DMPlexGetChart(dm, &pStart, &pEnd);
848:   for (p = pStart; p < pEnd; p++) {
849:     PetscInt bval;

851:     DMLabelGetValue(boundary, p, &bval);
852:     if (bval == 1) {
853:       PetscInt dof;

855:       PetscSectionGetDof(section, p, &dof);
856:       PetscSectionSetConstraintDof(section, p, dof);
857:     }
858:   }
859:   DMLabelDestroy(&boundary);
860:   PetscSectionSetUp(section);
861:   return 0;
862: }

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

867:   Collective on sp

869:   Input Parameters:
870: . sp      - The PetscDualSpace

872:   Output Parameter:
873: . section - The section

875:   Level: advanced

877: .seealso: PetscDualSpaceCreate(), DMPLEX
878: @*/
879: PetscErrorCode PetscDualSpaceGetSection(PetscDualSpace sp, PetscSection *section)
880: {
881:   PetscInt       pStart, pEnd, p;

883:   if (!sp->pointSection) {
884:     /* mark the boundary */
885:     PetscDualSpaceSectionCreate_Internal(sp, &(sp->pointSection));
886:     DMPlexGetChart(sp->dm,&pStart,&pEnd);
887:     for (p = pStart; p < pEnd; p++) {
888:       PetscDualSpace psp;

890:       PetscDualSpaceGetPointSubspace(sp, p, &psp);
891:       if (psp) {
892:         PetscInt dof;

894:         PetscDualSpaceGetInteriorDimension(psp, &dof);
895:         PetscSectionSetDof(sp->pointSection,p,dof);
896:       }
897:     }
898:     PetscDualSpaceSectionSetUp_Internal(sp,sp->pointSection);
899:   }
900:   *section = sp->pointSection;
901:   return 0;
902: }

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

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

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

939:       PetscDualSpaceGetFunctional(ssp, f, &fn);
940:       PetscQuadraturePushForward(fn, dim, sv0, v0, J, k, &(sp->functional[off+f]));
941:     }
942:   }
943:   PetscFree3(v0, sv0, J);
944:   return 0;
945: }

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

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

959:   Output Parameter:
960: . value   - numComp output values

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

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

967:   Level: beginner

969: .seealso: PetscDualSpaceCreate()
970: @*/
971: 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)
972: {
976:   (*sp->ops->apply)(sp, f, time, cgeom, numComp, func, ctx, value);
977:   return 0;
978: }

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

983:   Input Parameters:
984: + sp        - The PetscDualSpace object
985: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()

987:   Output Parameter:
988: . spValue   - The values of all dual space functionals

990:   Level: beginner

992: .seealso: PetscDualSpaceCreate()
993: @*/
994: PetscErrorCode PetscDualSpaceApplyAll(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
995: {
997:   (*sp->ops->applyall)(sp, pointEval, spValue);
998:   return 0;
999: }

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

1004:   Input Parameters:
1005: + sp        - The PetscDualSpace object
1006: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()

1008:   Output Parameter:
1009: . spValue   - The values of interior dual space functionals

1011:   Level: beginner

1013: .seealso: PetscDualSpaceCreate()
1014: @*/
1015: PetscErrorCode PetscDualSpaceApplyInterior(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1016: {
1018:   (*sp->ops->applyint)(sp, pointEval, spValue);
1019:   return 0;
1020: }

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

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

1034:   Output Parameter:
1035: . value   - The output value

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

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

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

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

1046: where both n and f have Nc components.

1048:   Level: beginner

1050: .seealso: PetscDualSpaceCreate()
1051: @*/
1052: 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)
1053: {
1054:   DM               dm;
1055:   PetscQuadrature  n;
1056:   const PetscReal *points, *weights;
1057:   PetscReal        x[3];
1058:   PetscScalar     *val;
1059:   PetscInt         dim, dE, qNc, c, Nq, q;
1060:   PetscBool        isAffine;

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

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

1091:   Input Parameters:
1092: + sp        - The PetscDualSpace object
1093: - pointEval - Evaluation at the points returned by PetscDualSpaceGetAllData()

1095:   Output Parameter:
1096: . spValue   - The values of all dual space functionals

1098:   Level: beginner

1100: .seealso: PetscDualSpaceCreate()
1101: @*/
1102: PetscErrorCode PetscDualSpaceApplyAllDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1103: {
1104:   Vec              pointValues, dofValues;
1105:   Mat              allMat;

1110:   PetscDualSpaceGetAllData(sp, NULL, &allMat);
1111:   if (!(sp->allNodeValues)) {
1112:     MatCreateVecs(allMat, &(sp->allNodeValues), NULL);
1113:   }
1114:   pointValues = sp->allNodeValues;
1115:   if (!(sp->allDofValues)) {
1116:     MatCreateVecs(allMat, NULL, &(sp->allDofValues));
1117:   }
1118:   dofValues = sp->allDofValues;
1119:   VecPlaceArray(pointValues, pointEval);
1120:   VecPlaceArray(dofValues, spValue);
1121:   MatMult(allMat, pointValues, dofValues);
1122:   VecResetArray(dofValues);
1123:   VecResetArray(pointValues);
1124:   return 0;
1125: }

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

1130:   Input Parameters:
1131: + sp        - The PetscDualSpace object
1132: - pointEval - Evaluation at the points returned by PetscDualSpaceGetInteriorData()

1134:   Output Parameter:
1135: . spValue   - The values of interior dual space functionals

1137:   Level: beginner

1139: .seealso: PetscDualSpaceCreate()
1140: @*/
1141: PetscErrorCode PetscDualSpaceApplyInteriorDefault(PetscDualSpace sp, const PetscScalar *pointEval, PetscScalar *spValue)
1142: {
1143:   Vec              pointValues, dofValues;
1144:   Mat              intMat;

1149:   PetscDualSpaceGetInteriorData(sp, NULL, &intMat);
1150:   if (!(sp->intNodeValues)) {
1151:     MatCreateVecs(intMat, &(sp->intNodeValues), NULL);
1152:   }
1153:   pointValues = sp->intNodeValues;
1154:   if (!(sp->intDofValues)) {
1155:     MatCreateVecs(intMat, NULL, &(sp->intDofValues));
1156:   }
1157:   dofValues = sp->intDofValues;
1158:   VecPlaceArray(pointValues, pointEval);
1159:   VecPlaceArray(dofValues, spValue);
1160:   MatMult(intMat, pointValues, dofValues);
1161:   VecResetArray(dofValues);
1162:   VecResetArray(pointValues);
1163:   return 0;
1164: }

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

1169:   Input Parameter:
1170: . sp - The dualspace

1172:   Output Parameters:
1173: + allNodes - A PetscQuadrature object containing all evaluation nodes
1174: - allMat - A Mat for the node-to-dof transformation

1176:   Level: advanced

1178: .seealso: PetscDualSpaceCreate()
1179: @*/
1180: PetscErrorCode PetscDualSpaceGetAllData(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1181: {
1185:   if ((!sp->allNodes || !sp->allMat) && sp->ops->createalldata) {
1186:     PetscQuadrature qpoints;
1187:     Mat amat;

1189:     (*sp->ops->createalldata)(sp,&qpoints,&amat);
1190:     PetscQuadratureDestroy(&(sp->allNodes));
1191:     MatDestroy(&(sp->allMat));
1192:     sp->allNodes = qpoints;
1193:     sp->allMat = amat;
1194:   }
1195:   if (allNodes) *allNodes = sp->allNodes;
1196:   if (allMat) *allMat = sp->allMat;
1197:   return 0;
1198: }

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

1203:   Input Parameter:
1204: . sp - The dualspace

1206:   Output Parameters:
1207: + allNodes - A PetscQuadrature object containing all evaluation nodes
1208: - allMat - A Mat for the node-to-dof transformation

1210:   Level: advanced

1212: .seealso: PetscDualSpaceCreate()
1213: @*/
1214: PetscErrorCode PetscDualSpaceCreateAllDataDefault(PetscDualSpace sp, PetscQuadrature *allNodes, Mat *allMat)
1215: {
1216:   PetscInt        spdim;
1217:   PetscInt        numPoints, offset;
1218:   PetscReal       *points;
1219:   PetscInt        f, dim;
1220:   PetscInt        Nc, nrows, ncols;
1221:   PetscInt        maxNumPoints;
1222:   PetscQuadrature q;
1223:   Mat             A;

1225:   PetscDualSpaceGetNumComponents(sp, &Nc);
1226:   PetscDualSpaceGetDimension(sp,&spdim);
1227:   if (!spdim) {
1228:     PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1229:     PetscQuadratureSetData(*allNodes,0,0,0,NULL,NULL);
1230:   }
1231:   nrows = spdim;
1232:   PetscDualSpaceGetFunctional(sp,0,&q);
1233:   PetscQuadratureGetData(q,&dim,NULL,&numPoints,NULL,NULL);
1234:   maxNumPoints = numPoints;
1235:   for (f = 1; f < spdim; f++) {
1236:     PetscInt Np;

1238:     PetscDualSpaceGetFunctional(sp,f,&q);
1239:     PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1240:     numPoints += Np;
1241:     maxNumPoints = PetscMax(maxNumPoints,Np);
1242:   }
1243:   ncols = numPoints * Nc;
1244:   PetscMalloc1(dim*numPoints,&points);
1245:   MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, maxNumPoints * Nc, NULL, &A);
1246:   for (f = 0, offset = 0; f < spdim; f++) {
1247:     const PetscReal *p, *w;
1248:     PetscInt        Np, i;
1249:     PetscInt        fnc;

1251:     PetscDualSpaceGetFunctional(sp,f,&q);
1252:     PetscQuadratureGetData(q,NULL,&fnc,&Np,&p,&w);
1254:     for (i = 0; i < Np * dim; i++) {
1255:       points[offset* dim + i] = p[i];
1256:     }
1257:     for (i = 0; i < Np * Nc; i++) {
1258:       MatSetValue(A, f, offset * Nc, w[i], INSERT_VALUES);
1259:     }
1260:     offset += Np;
1261:   }
1262:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1263:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1264:   PetscQuadratureCreate(PETSC_COMM_SELF,allNodes);
1265:   PetscQuadratureSetData(*allNodes,dim,0,numPoints,points,NULL);
1266:   *allMat = A;
1267:   return 0;
1268: }

1270: /*@
1271:   PetscDualSpaceGetInteriorData - Get all quadrature points necessary to compute the interior degrees of freedom from
1272:   this space, as well as the matrix that computes the degrees of freedom from the quadrature values.  Degrees of
1273:   freedom are interior degrees of freedom if they belong (by PetscDualSpaceGetSection()) to interior points in the
1274:   reference DMPlex: complementary boundary degrees of freedom are marked as constrained in the section returned by
1275:   PetscDualSpaceGetSection()).

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

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

1286:   Level: advanced

1288: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetDimension(), PetscDualSpaceGetNumComponents(), PetscQuadratureGetData()
1289: @*/
1290: PetscErrorCode PetscDualSpaceGetInteriorData(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1291: {
1295:   if ((!sp->intNodes || !sp->intMat) && sp->ops->createintdata) {
1296:     PetscQuadrature qpoints;
1297:     Mat imat;

1299:     (*sp->ops->createintdata)(sp,&qpoints,&imat);
1300:     PetscQuadratureDestroy(&(sp->intNodes));
1301:     MatDestroy(&(sp->intMat));
1302:     sp->intNodes = qpoints;
1303:     sp->intMat = imat;
1304:   }
1305:   if (intNodes) *intNodes = sp->intNodes;
1306:   if (intMat) *intMat = sp->intMat;
1307:   return 0;
1308: }

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

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

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

1322:   Level: advanced

1324: .seealso: PetscDualSpaceCreate(), PetscDualSpaceGetInteriorData()
1325: @*/
1326: PetscErrorCode PetscDualSpaceCreateInteriorDataDefault(PetscDualSpace sp, PetscQuadrature *intNodes, Mat *intMat)
1327: {
1328:   DM              dm;
1329:   PetscInt        spdim0;
1330:   PetscInt        Nc;
1331:   PetscInt        pStart, pEnd, p, f;
1332:   PetscSection    section;
1333:   PetscInt        numPoints, offset, matoffset;
1334:   PetscReal       *points;
1335:   PetscInt        dim;
1336:   PetscInt        *nnz;
1337:   PetscQuadrature q;
1338:   Mat             imat;

1341:   PetscDualSpaceGetSection(sp, &section);
1342:   PetscSectionGetConstrainedStorageSize(section, &spdim0);
1343:   if (!spdim0) {
1344:     *intNodes = NULL;
1345:     *intMat = NULL;
1346:     return 0;
1347:   }
1348:   PetscDualSpaceGetNumComponents(sp, &Nc);
1349:   PetscSectionGetChart(section, &pStart, &pEnd);
1350:   PetscDualSpaceGetDM(sp, &dm);
1351:   DMGetDimension(dm, &dim);
1352:   PetscMalloc1(spdim0, &nnz);
1353:   for (p = pStart, f = 0, numPoints = 0; p < pEnd; p++) {
1354:     PetscInt dof, cdof, off, d;

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

1363:       PetscDualSpaceGetFunctional(sp,off,&q);
1364:       PetscQuadratureGetData(q,NULL,NULL,&Np,NULL,NULL);
1365:       nnz[f] = Np * Nc;
1366:       numPoints += Np;
1367:     }
1368:   }
1369:   MatCreateSeqAIJ(PETSC_COMM_SELF, spdim0, numPoints * Nc, 0, nnz, &imat);
1370:   PetscFree(nnz);
1371:   PetscMalloc1(dim*numPoints,&points);
1372:   for (p = pStart, f = 0, offset = 0, matoffset = 0; p < pEnd; p++) {
1373:     PetscInt dof, cdof, off, d;

1375:     PetscSectionGetDof(section, p, &dof);
1376:     PetscSectionGetConstraintDof(section, p, &cdof);
1377:     if (!(dof - cdof)) continue;
1378:     PetscSectionGetOffset(section, p, &off);
1379:     for (d = 0; d < dof; d++, off++, f++) {
1380:       const PetscReal *p;
1381:       const PetscReal *w;
1382:       PetscInt        Np, i;

1384:       PetscDualSpaceGetFunctional(sp,off,&q);
1385:       PetscQuadratureGetData(q,NULL,NULL,&Np,&p,&w);
1386:       for (i = 0; i < Np * dim; i++) {
1387:         points[offset + i] = p[i];
1388:       }
1389:       for (i = 0; i < Np * Nc; i++) {
1390:         MatSetValue(imat, f, matoffset + i, w[i],INSERT_VALUES);
1391:       }
1392:       offset += Np * dim;
1393:       matoffset += Np * Nc;
1394:     }
1395:   }
1396:   PetscQuadratureCreate(PETSC_COMM_SELF,intNodes);
1397:   PetscQuadratureSetData(*intNodes,dim,0,numPoints,points,NULL);
1398:   MatAssemblyBegin(imat, MAT_FINAL_ASSEMBLY);
1399:   MatAssemblyEnd(imat, MAT_FINAL_ASSEMBLY);
1400:   *intMat = imat;
1401:   return 0;
1402: }

1404: /*@
1405:   PetscDualSpaceEqual - Determine if a dual space is equivalent

1407:   Input Parameters:
1408: + A    - A PetscDualSpace object
1409: - B    - Another PetscDualSpace object

1411:   Output Parameter:
1412: . equal - PETSC_TRUE if the dual spaces are equivalent

1414:   Level: advanced

1416: .seealso: PetscDualSpaceCreate()
1417: @*/
1418: PetscErrorCode PetscDualSpaceEqual(PetscDualSpace A, PetscDualSpace B, PetscBool *equal)
1419: {
1420:   PetscInt sizeA, sizeB, dimA, dimB;
1421:   const PetscInt *dofA, *dofB;
1422:   PetscQuadrature quadA, quadB;
1423:   Mat matA, matB;

1428:   *equal = PETSC_FALSE;
1429:   PetscDualSpaceGetDimension(A, &sizeA);
1430:   PetscDualSpaceGetDimension(B, &sizeB);
1431:   if (sizeB != sizeA) {
1432:     return 0;
1433:   }
1434:   DMGetDimension(A->dm, &dimA);
1435:   DMGetDimension(B->dm, &dimB);
1436:   if (dimA != dimB) {
1437:     return 0;
1438:   }

1440:   PetscDualSpaceGetNumDof(A, &dofA);
1441:   PetscDualSpaceGetNumDof(B, &dofB);
1442:   for (PetscInt d=0; d<dimA; d++) {
1443:     if (dofA[d] != dofB[d]) {
1444:       return 0;
1445:     }
1446:   }

1448:   PetscDualSpaceGetInteriorData(A, &quadA, &matA);
1449:   PetscDualSpaceGetInteriorData(B, &quadB, &matB);
1450:   if (!quadA && !quadB) {
1451:     *equal = PETSC_TRUE;
1452:   } else if (quadA && quadB) {
1453:     PetscQuadratureEqual(quadA, quadB, equal);
1454:     if (*equal == PETSC_FALSE) return 0;
1455:     if (!matA && !matB) return 0;
1456:     if (matA && matB) MatEqual(matA, matB, equal);
1457:     else *equal = PETSC_FALSE;
1458:   }
1459:   return 0;
1460: }

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

1465:   Input Parameters:
1466: + sp    - The PetscDualSpace object
1467: . f     - The basis functional index
1468: . time  - The time
1469: . cgeom - A context with geometric information for this cell, we currently just use the centroid
1470: . Nc    - The number of components for the function
1471: . func  - The input function
1472: - ctx   - A context for the function

1474:   Output Parameter:
1475: . value - The output value (scalar)

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

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

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

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

1486: where both n and f have Nc components.

1488:   Level: beginner

1490: .seealso: PetscDualSpaceCreate()
1491: @*/
1492: 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)
1493: {
1494:   DM               dm;
1495:   PetscQuadrature  n;
1496:   const PetscReal *points, *weights;
1497:   PetscScalar     *val;
1498:   PetscInt         dimEmbed, qNc, c, Nq, q;

1502:   PetscDualSpaceGetDM(sp, &dm);
1503:   DMGetCoordinateDim(dm, &dimEmbed);
1504:   PetscDualSpaceGetFunctional(sp, f, &n);
1505:   PetscQuadratureGetData(n, NULL, &qNc, &Nq, &points, &weights);
1507:   DMGetWorkArray(dm, Nc, MPIU_SCALAR, &val);
1508:   *value = 0.;
1509:   for (q = 0; q < Nq; ++q) {
1510:     (*func)(dimEmbed, time, cgeom->centroid, Nc, val, ctx);
1511:     for (c = 0; c < Nc; ++c) {
1512:       *value += val[c]*weights[q*Nc+c];
1513:     }
1514:   }
1515:   DMRestoreWorkArray(dm, Nc, MPIU_SCALAR, &val);
1516:   return 0;
1517: }

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

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

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

1529:   Not collective

1531:   Input Parameters:
1532: + sp - the PetscDualSpace object
1533: - height - the height of the mesh point for which the subspace is desired

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

1539:   Level: advanced

1541: .seealso: PetscSpaceGetHeightSubspace(), PetscDualSpace
1542: @*/
1543: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *subsp)
1544: {
1545:   PetscInt       depth = -1, cStart, cEnd;
1546:   DM             dm;

1551:   *subsp = NULL;
1552:   dm = sp->dm;
1553:   DMPlexGetDepth(dm, &depth);
1555:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1556:   if (height == 0 && cEnd == cStart + 1) {
1557:     *subsp = sp;
1558:     return 0;
1559:   }
1560:   if (!sp->heightSpaces) {
1561:     PetscInt h;
1562:     PetscCalloc1(depth+1, &(sp->heightSpaces));

1564:     for (h = 0; h <= depth; h++) {
1565:       if (h == 0 && cEnd == cStart + 1) continue;
1566:       if (sp->ops->createheightsubspace) (*sp->ops->createheightsubspace)(sp,height,&(sp->heightSpaces[h]));
1567:       else if (sp->pointSpaces) {
1568:         PetscInt hStart, hEnd;

1570:         DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
1571:         if (hEnd > hStart) {
1572:           const char *name;

1574:           PetscObjectReference((PetscObject)(sp->pointSpaces[hStart]));
1575:           if (sp->pointSpaces[hStart]) {
1576:             PetscObjectGetName((PetscObject) sp,                     &name);
1577:             PetscObjectSetName((PetscObject) sp->pointSpaces[hStart], name);
1578:           }
1579:           sp->heightSpaces[h] = sp->pointSpaces[hStart];
1580:         }
1581:       }
1582:     }
1583:   }
1584:   *subsp = sp->heightSpaces[height];
1585:   return 0;
1586: }

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

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

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

1597:   Not collective

1599:   Input Parameters:
1600: + sp - the PetscDualSpace object
1601: - point - the point (in the dual space's DM) for which the subspace is desired

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

1607:   Level: advanced

1609: .seealso: PetscDualSpace
1610: @*/
1611: PetscErrorCode PetscDualSpaceGetPointSubspace(PetscDualSpace sp, PetscInt point, PetscDualSpace *bdsp)
1612: {
1613:   PetscInt       pStart = 0, pEnd = 0, cStart, cEnd;
1614:   DM             dm;

1618:   *bdsp = NULL;
1619:   dm = sp->dm;
1620:   DMPlexGetChart(dm, &pStart, &pEnd);
1622:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1623:   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 */
1624:     *bdsp = sp;
1625:     return 0;
1626:   }
1627:   if (!sp->pointSpaces) {
1628:     PetscInt p;
1629:     PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));

1631:     for (p = 0; p < pEnd - pStart; p++) {
1632:       if (p + pStart == cStart && cEnd == cStart + 1) continue;
1633:       if (sp->ops->createpointsubspace) (*sp->ops->createpointsubspace)(sp,p+pStart,&(sp->pointSpaces[p]));
1634:       else if (sp->heightSpaces || sp->ops->createheightsubspace) {
1635:         PetscInt dim, depth, height;
1636:         DMLabel  label;

1638:         DMPlexGetDepth(dm,&dim);
1639:         DMPlexGetDepthLabel(dm,&label);
1640:         DMLabelGetValue(label,p+pStart,&depth);
1641:         height = dim - depth;
1642:         PetscDualSpaceGetHeightSubspace(sp, height, &(sp->pointSpaces[p]));
1643:         PetscObjectReference((PetscObject)sp->pointSpaces[p]);
1644:       }
1645:     }
1646:   }
1647:   *bdsp = sp->pointSpaces[point - pStart];
1648:   return 0;
1649: }

1651: /*@C
1652:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

1654:   Not collective

1656:   Input Parameter:
1657: . sp - the PetscDualSpace object

1659:   Output Parameters:
1660: + perms - Permutations of the interior degrees of freedom, parameterized by the point orientation
1661: - flips - Sign reversal of the interior degrees of freedom, parameterized by the point orientation

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

1667:   Level: developer

1669: @*/
1670: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
1671: {
1675:   if (sp->ops->getsymmetries) (sp->ops->getsymmetries)(sp,perms,flips);
1676:   return 0;
1677: }

1679: /*@
1680:   PetscDualSpaceGetFormDegree - Get the form degree k for the k-form the describes the pushforwards/pullbacks of this
1681:   dual space's functionals.

1683:   Input Parameter:
1684: . dsp - The PetscDualSpace

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

1694:   Level: developer

1696: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1697: @*/
1698: PetscErrorCode PetscDualSpaceGetFormDegree(PetscDualSpace dsp, PetscInt *k)
1699: {
1703:   *k = dsp->k;
1704:   return 0;
1705: }

1707: /*@
1708:   PetscDualSpaceSetFormDegree - Set the form degree k for the k-form the describes the pushforwards/pullbacks of this
1709:   dual space's functionals.

1711:   Input Parameters:
1712: + dsp - The PetscDualSpace
1713: - k   - The *signed* degree k of the k.  If k >= 0, this means that the degrees of freedom are k-forms, and are stored
1714:         in lexicographic order according to the basis of k-forms constructed from the wedge product of 1-forms.  So for example,
1715:         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).
1716:         If k < 0, this means that the degrees transform as k-forms, but are stored as (N-k) forms according to the
1717:         Hodge star map.  So for example if k = -2 and N = 3, this means that the degrees of freedom transform as 2-forms
1718:         but are stored as 1-forms.

1720:   Level: developer

1722: .seealso: PetscDTAltV, PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1723: @*/
1724: PetscErrorCode PetscDualSpaceSetFormDegree(PetscDualSpace dsp, PetscInt k)
1725: {
1726:   PetscInt dim;

1731:   dim = dsp->dm->dim;
1733:   dsp->k = k;
1734:   return 0;
1735: }

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

1740:   Input Parameter:
1741: . dsp - The PetscDualSpace

1743:   Output Parameter:
1744: . k   - The simplex dimension

1746:   Level: developer

1748:   Note: Currently supported values are
1749: $ 0: These are H_1 methods that only transform coordinates
1750: $ 1: These are Hcurl methods that transform functions using the covariant Piola transform (COVARIANT_PIOLA_TRANSFORM)
1751: $ 2: These are the same as 1
1752: $ 3: These are Hdiv methods that transform functions using the contravariant Piola transform (CONTRAVARIANT_PIOLA_TRANSFORM)

1754: .seealso: PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceTransformType
1755: @*/
1756: PetscErrorCode PetscDualSpaceGetDeRahm(PetscDualSpace dsp, PetscInt *k)
1757: {
1758:   PetscInt dim;

1763:   dim = dsp->dm->dim;
1764:   if (!dsp->k) *k = IDENTITY_TRANSFORM;
1765:   else if (dsp->k == 1) *k = COVARIANT_PIOLA_TRANSFORM;
1766:   else if (dsp->k == -(dim - 1)) *k = CONTRAVARIANT_PIOLA_TRANSFORM;
1767:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported transformation");
1768:   return 0;
1769: }

1771: /*@C
1772:   PetscDualSpaceTransform - Transform the function values

1774:   Input Parameters:
1775: + dsp       - The PetscDualSpace
1776: . trans     - The type of transform
1777: . isInverse - Flag to invert the transform
1778: . fegeom    - The cell geometry
1779: . Nv        - The number of function samples
1780: . Nc        - The number of function components
1781: - vals      - The function values

1783:   Output Parameter:
1784: . vals      - The transformed function values

1786:   Level: intermediate

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

1790: .seealso: PetscDualSpaceTransformGradient(), PetscDualSpaceTransformHessian(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1791: @*/
1792: PetscErrorCode PetscDualSpaceTransform(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1793: {
1794:   PetscReal Jstar[9] = {0};
1795:   PetscInt dim, v, c, Nk;

1801:   /* TODO: not handling dimEmbed != dim right now */
1802:   dim = dsp->dm->dim;
1803:   /* No change needed for 0-forms */
1804:   if (!dsp->k) return 0;
1805:   PetscDTBinomialInt(dim, PetscAbsInt(dsp->k), &Nk);
1806:   /* TODO: use fegeom->isAffine */
1807:   PetscDTAltVPullbackMatrix(dim, dim, isInverse ? fegeom->J : fegeom->invJ, dsp->k, Jstar);
1808:   for (v = 0; v < Nv; ++v) {
1809:     switch (Nk) {
1810:     case 1:
1811:       for (c = 0; c < Nc; c++) vals[v*Nc + c] *= Jstar[0];
1812:       break;
1813:     case 2:
1814:       for (c = 0; c < Nc; c += 2) DMPlex_Mult2DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1815:       break;
1816:     case 3:
1817:       for (c = 0; c < Nc; c += 3) DMPlex_Mult3DReal_Internal(Jstar, 1, &vals[v*Nc + c], &vals[v*Nc + c]);
1818:       break;
1819:     default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported form size %D for transformation", Nk);
1820:     }
1821:   }
1822:   return 0;
1823: }

1825: /*@C
1826:   PetscDualSpaceTransformGradient - Transform the function gradient values

1828:   Input Parameters:
1829: + dsp       - The PetscDualSpace
1830: . trans     - The type of transform
1831: . isInverse - Flag to invert the transform
1832: . fegeom    - The cell geometry
1833: . Nv        - The number of function gradient samples
1834: . Nc        - The number of function components
1835: - vals      - The function gradient values

1837:   Output Parameter:
1838: . vals      - The transformed function gradient values

1840:   Level: intermediate

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

1844: .seealso: PetscDualSpaceTransform(), PetscDualSpacePullback(), PetscDualSpacePushforward(), PetscDualSpaceTransformType
1845: @*/
1846: PetscErrorCode PetscDualSpaceTransformGradient(PetscDualSpace dsp, PetscDualSpaceTransformType trans, PetscBool isInverse, PetscFEGeom *fegeom, PetscInt Nv, PetscInt Nc, PetscScalar vals[])
1847: {
1848:   const PetscInt dim = dsp->dm->dim, dE = fegeom->dimEmbed;
1849:   PetscInt       v, c, d;

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

1938: /*@C
1939:   PetscDualSpaceTransformHessian - Transform the function Hessian values

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

1950:   Output Parameter:
1951: . vals      - The transformed function Hessian values

1953:   Level: intermediate

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

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

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

2003: /*@C
2004:   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.

2006:   Input Parameters:
2007: + dsp        - The PetscDualSpace
2008: . fegeom     - The geometry for this cell
2009: . Nq         - The number of function samples
2010: . Nc         - The number of function components
2011: - pointEval  - The function values

2013:   Output Parameter:
2014: . pointEval  - The transformed function values

2016:   Level: advanced

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

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

2022: .seealso: PetscDualSpacePushforward(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2023: @*/
2024: PetscErrorCode PetscDualSpacePullback(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2025: {
2026:   PetscDualSpaceTransformType trans;
2027:   PetscInt                    k;

2033:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2034:      This determines their transformation properties. */
2035:   PetscDualSpaceGetDeRahm(dsp, &k);
2036:   switch (k)
2037:   {
2038:     case 0: /* H^1 point evaluations */
2039:     trans = IDENTITY_TRANSFORM;break;
2040:     case 1: /* Hcurl preserves tangential edge traces  */
2041:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2042:     case 2:
2043:     case 3: /* Hdiv preserve normal traces */
2044:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2045:     default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2046:   }
2047:   PetscDualSpaceTransform(dsp, trans, PETSC_TRUE, fegeom, Nq, Nc, pointEval);
2048:   return 0;
2049: }

2051: /*@C
2052:   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.

2054:   Input Parameters:
2055: + dsp        - The PetscDualSpace
2056: . fegeom     - The geometry for this cell
2057: . Nq         - The number of function samples
2058: . Nc         - The number of function components
2059: - pointEval  - The function values

2061:   Output Parameter:
2062: . pointEval  - The transformed function values

2064:   Level: advanced

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

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

2070: .seealso: PetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2071: @*/
2072: PetscErrorCode PetscDualSpacePushforward(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2073: {
2074:   PetscDualSpaceTransformType trans;
2075:   PetscInt                    k;

2081:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2082:      This determines their transformation properties. */
2083:   PetscDualSpaceGetDeRahm(dsp, &k);
2084:   switch (k)
2085:   {
2086:     case 0: /* H^1 point evaluations */
2087:     trans = IDENTITY_TRANSFORM;break;
2088:     case 1: /* Hcurl preserves tangential edge traces  */
2089:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2090:     case 2:
2091:     case 3: /* Hdiv preserve normal traces */
2092:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2093:     default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2094:   }
2095:   PetscDualSpaceTransform(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2096:   return 0;
2097: }

2099: /*@C
2100:   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.

2102:   Input Parameters:
2103: + dsp        - The PetscDualSpace
2104: . fegeom     - The geometry for this cell
2105: . Nq         - The number of function gradient samples
2106: . Nc         - The number of function components
2107: - pointEval  - The function gradient values

2109:   Output Parameter:
2110: . pointEval  - The transformed function gradient values

2112:   Level: advanced

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

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

2118: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2119: @*/
2120: PetscErrorCode PetscDualSpacePushforwardGradient(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2121: {
2122:   PetscDualSpaceTransformType trans;
2123:   PetscInt                    k;

2129:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2130:      This determines their transformation properties. */
2131:   PetscDualSpaceGetDeRahm(dsp, &k);
2132:   switch (k)
2133:   {
2134:     case 0: /* H^1 point evaluations */
2135:     trans = IDENTITY_TRANSFORM;break;
2136:     case 1: /* Hcurl preserves tangential edge traces  */
2137:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2138:     case 2:
2139:     case 3: /* Hdiv preserve normal traces */
2140:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2141:     default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2142:   }
2143:   PetscDualSpaceTransformGradient(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2144:   return 0;
2145: }

2147: /*@C
2148:   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.

2150:   Input Parameters:
2151: + dsp        - The PetscDualSpace
2152: . fegeom     - The geometry for this cell
2153: . Nq         - The number of function Hessian samples
2154: . Nc         - The number of function components
2155: - pointEval  - The function gradient values

2157:   Output Parameter:
2158: . pointEval  - The transformed function Hessian values

2160:   Level: advanced

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

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

2166: .seealso: PetscDualSpacePushforward(), PPetscDualSpacePullback(), PetscDualSpaceTransform(), PetscDualSpaceGetDeRahm()
2167: @*/
2168: PetscErrorCode PetscDualSpacePushforwardHessian(PetscDualSpace dsp, PetscFEGeom *fegeom, PetscInt Nq, PetscInt Nc, PetscScalar pointEval[])
2169: {
2170:   PetscDualSpaceTransformType trans;
2171:   PetscInt                    k;

2177:   /* The dualspace dofs correspond to some simplex in the DeRahm complex, which we label by k.
2178:      This determines their transformation properties. */
2179:   PetscDualSpaceGetDeRahm(dsp, &k);
2180:   switch (k)
2181:   {
2182:     case 0: /* H^1 point evaluations */
2183:     trans = IDENTITY_TRANSFORM;break;
2184:     case 1: /* Hcurl preserves tangential edge traces  */
2185:     trans = COVARIANT_PIOLA_TRANSFORM;break;
2186:     case 2:
2187:     case 3: /* Hdiv preserve normal traces */
2188:     trans = CONTRAVARIANT_PIOLA_TRANSFORM;break;
2189:     default: SETERRQ(PetscObjectComm((PetscObject) dsp), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported simplex dim %D for transformation", k);
2190:   }
2191:   PetscDualSpaceTransformHessian(dsp, trans, PETSC_FALSE, fegeom, Nq, Nc, pointEval);
2192:   return 0;
2193: }